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Abstract 

System  identification  has  long  been  used  as  a  tool  for  flight  test  engineers  to  char¬ 
acterize  systems  under  test;  however,  the  inputs  to  these  characterization  activities 
have  previously  been  limited  to  wind  tunnel  and  flight  test  data.  There  has  been 
a  recent  effort  to  incorporate  computational  fluid  dynamics  (CFD)  into  the  system 
identification  process.  An  integral  piece  of  the  process  is  the  simulation  of  training 
maneuvers  utilizing  CFD.  Up  until  now,  the  suitability  of  particular  training  ma¬ 
neuvers  has  been  assessed  by  comparing  the  reduced  model  to  known  results.  This 
research  strives  to  recommend  a  set  of  parameters  to  be  used  in  determining  a  priori 
whether  a  training  maneuver  will  be  suitable  under  a  particular  set  of  flow  conditions. 
As  part  of  this  study,  regressor  space  parameters  (RSPs)  were  proposed  and  evalu¬ 
ated  for  several  maneuvers.  Then,  multivariate  polynomial  models  for  coefficient  of 
lift,  coefficient  of  drag  and  pitch  moment  coefficient  were  generated  from  the  results 
of  the  training  maneuvers  simulated  via  CREATE-AV’s  flow  solver,  Kestrel.  The 
resulting  model-predicted  coefficient  values  were  then  compared  to  CFD  simulations 
of  carefully  selected  comparison  maneuvers,  as  measured  by  a  set  of  suitable  statis¬ 
tical  metrics.  Using  the  statistical  analysis  tool  JMP®,  the  set  of  proposed  RSPs 
were  then  analyzed  in  order  to  evaluate  their  utility  in  assessing  the  suitability  of 
the  training  maneuvers.  This  research  showed  the  reliability  of  certain  RSPs  to  char¬ 
acterize  training  maneuver  performance.  Finally,  lessons  learned  from  this  research 
were  applied  to  improve  upon  current  best  practice  training  maneuvers. 
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TRAINING  MANEUVER  EVALUATION  FOR  REDUCED  ORDER  MODELING 


OF  STABILITY  &  CONTROL  PROPERTIES  USING  COMPUTATIONAL 

FLUID  DYNAMICS 

I.  Introduction 


1.1  General  Issue 

Design  changes  late  in  the  development  cycle  are  a  problem  that  have  plagued  mil¬ 
itary  fighter  aircraft  for  decades  [15].  While  analytical  and  numerical  techniques  can 
provide  some  insight  into  stability  and  control  (S&C)  characteristics  at  benign  flight 
conditions,  advanced  flight  conditions  that  include  non-linearities  are  not  easily  cal¬ 
culated.  Running  the  whole  spectrum  of  dynamic  maneuvers  using  a  fully-turbulent 
Navier-Stokes  solver  in  order  to  discover  possible  problems  in  the  S&C  characteristics 
is  too  costly. 

Specific  examples  of  costly  late  design  changes  are  numerous.  Consider  the  F- 
4,  at  high  a  (angle  of  attack);  the  aircraft  lost  directional  stability.  Over  200  F-4s 
were  lost  before  a  solution  was  found  [14],  The  737NG  had  undesirable  flaps-up  stall 
characteristics;  the  aircraft  experienced  stick  lightening  during  stall  [14],  Northrop 
Grumman  observed  yaw  departures  and  post  stall  lateral  and  directional  instability 
with  the  F-5  which  resulted  in  the  loss  of  a  flight  test  aircraft  and  large  cost  overruns 
and  schedule  delays  [14], 

There  is  a  Department  of  Defense  (DoD)  and  Air  Force  (AF)  need  to  find  and 
mitigate  these  risks  earlier  in  the  design  cycle  in  order  to  save  time,  money  and 
possibly  human  life. 
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1.2  Background 


There  are  several  avenues  for  obtaining  the  stability  characteristics  for  an  aircraft: 
analytical  calculations,  wind-tunnel  testing,  flight  testing,  and  full  computational 
simulations.  System  identification  SID  has  been  offered  as  methodology  for  utilizing 
computational  simulations  in  order  to  produce  simpler  numerical  calculations  via 
reduced  order  models  (ROMs)  to  calculate  S&C  characteristics.  The  paragraphs 
below  will  briefly  discuss  the  benefits,  challenges  and  limitations  for  each  method. 

Analytical  results  are  useful  in  that  the  engineer  can  begin  as  early  as  desired 
in  the  design  cycle.  To  some  degree,  analytical  results  guide  aircraft  designers  early 
in  the  design;  however,  there  is  a  limit  to  how  much  realism  analytical  techniques 
are  able  to  incorporate.  Non-linearities  and  complex  configurations  make  analytical 
techniques  an  impractical  choice. 

Wind  tunnel  testing  has  the  advantage  of  being  able  to  produce  a  wide  variety  of 
aerodynamic  results  as  soon  as  a  configuration  is  determined.  Unfortunately,  there 
are  also  limitations  to  the  results  wind  tunnel  testing  can  obtain.  A  wind  tunnel’s 
testing  parameters  are  unlikely  to  span  the  entire  operational  envelope  of  a  fighter 
aircraft.  Coefficients  such  as  coefficient  of  lift  (Cl)  and  coefficient  of  drag  (Co)  are 
easy  to  obtain  for  static  flight  conditions;  the  ability  to  evaluate  these  coefficients 
for  dynamic  maneuvers  is  limited.  This  difficultly  reduces  the  applicability  of  wind 
tunnel  testing  to  obtaining  derivatives  such  as  Cn&  which  is  the  partial  derivative  of 
Cjsr ,  to  cc. 

The  benefits  to  flight  testing  are  numerous.  While  the  data  provided  by  flight 
testing  are  absolutely  invaluable,  the  downside  is  that  a  fully  functional  prototype 
must  be  used,  and  the  prototype  will  likely  only  be  available  very  near  the  end  of 
the  design  process.  At  that  point,  considerable  time,  effort  and  money  have  been 
invested  in  the  project.  Major  design  changes  at  a  late  stage  are  detrimental  to 
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the  project,  particularly  on  modern  aircraft,  where  stealth  characteristics  are  highly 
dependent  on  aircraft  geometry.  If  sudden  instabilities  are  found  during  the  test,  it 
can  be  dangerous  for  the  aircraft  as  well  as  the  pilot.  Additionally,  flight  test  has  the 
downside  of  being  the  most  expensive  of  all  the  methods. 

Computational  results  also  have  the  benefit  of  yielding  results  once  a  geometry  is 
decided  upon-after  the  necessary  grid  generation,  grid  refinement  study  and  solution 
time.  However,  there  are  limitations  to  computational  fluid  dynamics  (CFD)  that 
need  to  be  understood  and  taken  into  account  in  order  to  gain  useful  results.  The 
quality  of  the  outputs  is  heavily  dependent  upon  the  quality  of  the  inputs.  Grid  and 
geometry  need  to  be  suitable;  otherwise,  no  matter  how  highly  accurate  the  solver 
may  be,  the  results  will  be  inaccurate  or  misleading.  Also,  while  CFD  has  come  a 
long  way  in  the  recent  past,  there  is  still  no  “perfect  solution.”  Any  computational 
simulation  is  going  to  include  a  variety  of  assumptions,  the  biggest  of  which  will  be 
the  model  used  to  simplify  the  turbulent  structures  in  the  flow.  One  benefit  is  that 
CFD  has  the  capability  to  provide  insight  into  aircraft  behavior  much  earlier  in  the 
design  cycle.  Finally,  there  is  also  the  benefit  of  having  data  everywhere  in  the  flow, 
as  opposed  to  a  few  select  areas  with  sensors  as  in  the  case  of  wind  tunnel  and  flight 
testing. 

Put  simply,  SID  is  a  process  by  which,  using  known  inputs  and  measured  outputs, 
a  complex  system  is  reduced  to  a  simplified  model.  In  this  work,  the  known  inputs  are 
parameters  such  as  a  and  Q.  The  measured  outputs  will  be  the  stability  coefficients: 
Cl,  Cd  and  Cm-  Previously,  system  identification  had  been  considered  a  useful 
process  only  for  flight  test  and  wind  tunnel  testing  data.  From  System  Identification, 
Theory  and  Practice :  “There  are  also  analytic  methods  to  find  the  aerodynamic 
function  dependencies,  such  as  computational  fluid  dynamics  (CFD). ..although  these 
methods  work  well  in  some  cases,  typically  for  low  angles  of  attack  and  low  rotational 
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rates,  the  best  aerodynamic  predictions  are  obtained  using  experimental  methods.” 
[20]  However,  there  has  been  a  recent  effort  to  incorporate  the  usefulness  of  the 
reduced  order  models  obtained  from  the  system  identification  process  utilizing  the 
early  availability  of  solutions  gained  from  high  hdelity-CFD. 

1.3  Research  Objectives 

This  project  seeks  to  further  the  current  understanding  of  using  CFD  in  the  system 
identification  process  in  order  to  gain  S&C  insight  earlier  in  the  design  process  and  to 
refine  the  methods  to  do  so.  In  order  for  CFD  to  be  used  in  the  system  identification 
process,  one  must  simulate  an  appropriate  training  maneuver.  First,  it  is  necessary 
to  determine  what  makes  a  good  training  maneuver.  Furthermore,  is  there  a  way  to 
evaluate  the  suitability  of  a  training  maneuver  a  priori ? 

The  focus  of  this  project  is  to  determine  a  set  of  parameters  that  will  indicate  if  a 
training  maneuver  is  suitable  prior  to  simulating  the  maneuver.  Previous  research  by 
Butler  [5],  has  been  conducted  in  order  to  find  parameters  to  be  used  prior  to  running 
a  maneuver.  However,  the  parameters  were  not  throughly  evaluated.  The  objective 
of  this  research  is  to  refine  and  create  parameters  in  order  to  estimate  future  model 
prediction  accuracy. 

The  path  ahead  for  this  project  consists  of  several  stages.  The  first  stage  is  to 
take  the  proposed  parameters  from  Butler  [5]  and  use  them  to  assess  the  suitability  of 
benchmark  maneuvers.  The  next  stage  is  to  simulate  the  different  training  maneuvers 
using  CFD  and  create  ROMs  from  the  outputs  using  SID  software.  The  next  step  is 
to  conduct  CFD  simulations  of  the  comparison  maneuvers  and  compare  the  model- 
predicted  coefficients  to  the  CFD  coefficients  using  selected  metrics.  Once  that  is 
completed,  the  next  step  is  to  evaluate  how  effectively  the  regressor  space  parameters 
(RSPs)  measured  the  prediction  capability  of  the  models.  If  necessary,  optimize 
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and  generate  new  parameters  to  improve  the  prediction  capability  of  the  models. 
Once  satisfied  with  the  parameters,  the  final  step  is  to  create  an  ideal  maneuver  that 
optimizes  those  RSPs. 

1.4  Road  Map 

Chapter  1  discussed  a  brief  overview  of  the  background  and  objectives  to  this  re¬ 
search.  Chapter  If  will  discuss  in  more  detail  the  background  information  for  this 
research.  Topics  will  include  the  governing  equations,  turbulence  theory,  turbulence 
modeling,  stability  &  control,  system  identification,  previous  research  and  Kestrel. 
Chapter  111  will  discuss  the  overall  methodology  used  in  the  research.  Topics  in 
Chapter  111  include  the  grid  convergence  study  method,  initial  regressor  space  pa¬ 
rameters,  training  maneuvers,  comparison  maneuvers,  reduced  order  model  genera¬ 
tion  and  analysis.  Chapter  IV  will  discuss  the  results  of  the  methods  employed  in 
Chapter  III,  to  include  the  grid  convergence  study,  the  regressor  space  parameters. 
Chapter  IV  will  detail  the  comparison  of  the  reduced  order  model  predictions  vs  the 
CFD  simulated  results,  and  how  those  results  motivated  the  design  of  two  additional 
training  maneuvers.  Chapter  V  will  tie  all  the  results  together  to  make  conclusions. 
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II.  Background 


This  chapter  will  present  a  discussion  of  necessary  background  material  for  this 
research.  Computational  fluid  dynamics  is  a  major  component  of  this  research,  and 
as  such,  the  background  includes  a  discussion  on  the  governing  equations  for  CFD  as 
well  as  turbulence  modeling  and  the  expected  flow  features  to  be  modeled.  Also,  this 
chapter  will  include  a  discussion  on  system  identification  and  the  process  by  which  a 
reduced  order  model  is  created.  Finally,  previous  research  in  this  area  of  study  will 
be  examined,  and  its  application  to  this  research  will  be  discussed. 

2.1  Computational  Fluid  Dynamics 

Computational  fluid  dynamics  is  the  combination  of  physics,  numerical  calcula¬ 
tions,  and  computer  science  with  the  goal  of  simulating  fluid  flows.  In  its  infancy, 
CFD  was  used  to  solve  the  2-D  Euler  equations.  With  the  increase  in  computing 
power  and  technology,  the  simulations  have  become  more  and  more  advanced.  Us¬ 
ing  multigrid,  multiblock,  implicit  time  integration  and  unstructured  grids,  CFD  has 
evolved  to  simulate  highly  unsteady,  viscous  flows  around  complex  geometries.  The 
following  discussion  will  focus  on  the  finite  volume  discretization  of  the  Navier-Stokes 
equations  with  turbulence  modeling,  as  that  is  the  focus  of  this  research. 

2.1.1  Governing  Equations 

There  are  three  governing  relations  for  general  CFD  computations,  which  are 
further  simplified  depending  on  the  fidelity  of  the  solution.  The  three  relations,  as 
shown  in  Blazek  [4],  are  the  conservation  of  mass,  conservation  of  momentum,  and 
conservation  of  energy  equations.  In  finite  volume  methods,  the  domain  is  discretized 
as  cells,  and  the  amount  of  the  conserved  quantities  crossing  the  boundaries  is  known 
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as  the  flux.  Taking  the  three  conservation  relations  in  integral  form  and  combining 
them  results  in  the  integral  form  of  the  Navier-Stokes  equations 


d 
d  t 


Fv)dS 


(2.1.1) 


where  W  is  a  vector  of  conserved  variables,  Fc  is  a  vector  of  convective  fluxes,  Fv  is  a 
vector  of  viscous  fluxes,  Q  is  a  collection  of  source  terms,  dFl  is  the  differential  volume 
and  dS  is  the  differential  surface. 

The  conservative  variables  are  given  in  the  vector 
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where  p  is  the  flow  density,  u  is  the  x-direction  velocity,  v  is  the  y-direction  velocity,  w 
is  the  z-direction  velocity  and  E  is  the  energy.  The  convective  fluxes  are  represented 
by  the  vector 
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where  V  —  U  ■  n  and  is  the  contravariant  velocity.  The  viscous  fluxes  are  given  by 
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where  is  the  sheer  stress  tensor.  Lastly,  the  source  terms  are  given  by 
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2.1.2  Turbulence  Theory  and  Modeling 

A  rigorous  definition  of  turbulence  is  hard  to  give;  however,  there  are  some  distinct 
characteristics  to  help  identify  turbulent  flow.  Unsteady,  irregular,  and  apparently 
random  fluctuations  are  three  distinct  characteristics.  Flow  is  likely  to  make  the 
switch  to  turbulence  when  small  perturbations  are  amplified.  The  sources  of  these 


small  perturbations  can  be  free  stream  turbulence,  surface  roughness  and  vibrations. 
The  transition  to  turbulence  depends  on  the  Reynolds  number.  At  low  Reynolds 
number,  the  viscous  forces  damp  out  the  instabilities  caused  by  turbulence.  At  high 
Reynolds  numbers,  inertial  disturbances  are  sufficiently  damped  and  none  are  ampli¬ 
fied.  Another  characteristic  of  turbulent  flow  is  the  large  range  of  length  and  time 
scales.  This  is  due  to  the  energy  cascade  in  turbulence.  At  the  largest  scales,  the 
inertial  effects  dominate  the  flow  behavior,  but  then  large  scale  vortices  are  stretched 
by  the  flow  gradients  which  transfers  energy  to  smaller  scales.  When  the  smallest 
scale  is  reached,  the  energy  is  dissipated  as  heat  due  to  viscosity-dominated  eddies. 
Kolmogorov  assumed  that  the  smallest  scales  are  universal  throughout  different  flows 
and  that  the  equilibrium  between  energy  states  from  the  large  turbulence  scales  can 
be  modeled  as 

dh 

(2.1.7) 


e  = 


dt  ’ 


where  e  is  the  TKE  dissipative  rate  and  k  is  the  turbulent  kinetic  energy.  The 
Kolmogorov  length  scale  is  expressed  as 


r3V/4 

Vk  =1  —  1 


e  J 


(2.1.8) 


Close  to  the  wall,  the  viscous  effects  dominate  and  u+  =  y+,  this  relation  is  the  law 
of  the  wall.  While  y+  <  5  that  region  is  known  as  the  viscous  sublayer.  Then  there 
is  the  log  layer  which  is  shown  by  equation  (2.1.9). 

The  log  law  states  that  the  average  velocity  of  the  flow  is  proportional  to  the  log 
of  the  distance  from  the  wall.  The  log  wall  is  given  by 


u+  =  -log(y+)  +  /3. 


(2.1.9) 
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where  k  is  the  Karman  constant  and  f3  is  another  constant.  The  various  sublayers 
are  shown  in  Figure  1. 


y+ 
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Figure  1.  Mean  velocity  distribution  across  a  turbulent  boundary  layer  [6] 

There  are  several  different  ways  to  model  turbulence.  Within  these  turbulence 
models  there  are  several  levels  of  approximation.  The  highest  fidelity  is  direct  numer¬ 
ical  simulation  (DNS).  In  decreasing  order  of  fidelity,  there  are  large  eddy  simulation 
(LES)  methods,  and  there  are  Reynolds  averaged  Navier-Stokes  (RANS)  methods. 
Bridging  the  gap  between  RANS  and  LES  are  hybrid  methods,  the  most  popular 
being  delayed  detached  eddy  simulation  (DDES).  The  range  of  turbulence  modeling 
techniques  is  shown  in  Figure  2 

DNS  is  the  most  accurate  of  all  types  of  turbulence  models,  as  it  resolves  all  of 
the  turbulent  length  scales.  This  direct  simulation  of  the  smallest  turbulence  scales 
requires  that  the  grid  must  be  able  to  resolve  down  to  the  Kolmogorov  length  scale 
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Figure  2.  Turbulence  modeling  approaches  [9] 

shown  in  Equation  (2.1.8)  [4],  The  small  grid  spacing  that  results  translates  to  an 
extremely  large  computational  cost,  which  causes  DNS  to  be  unpractical  for  most 
uses.  [9] 

In  LES,  the  flow  is  split  into  two  categories,  small  and  large  length  scales.  LES 
will  model  the  small  scales,  or  subgrid-scales,  while  simulating  the  large  length  scales. 
Since  LES  does  not  simulate  the  small  length  scales,  it  does  not  have  the  same  com¬ 
putational  costs  as  DNS.  However,  in  order  to  accurately  resolve  the  boundary  layer, 
the  grid  spacing  must  be  small.  The  grid  scales  must  be  so  small  that  it  currently 
drives  the  computational  costs  of  LES  solutions  out  of  general  engineering  use.  [9] 
RANS  models  substitute  an  average  and  a  fluctuating  value  into  the  N-S  equations 
and  then  make  some  qualifying  assumptions.  An  example  of  the  decomposed  values 
for  incompressible  flow  is 

Ui  =  Ui  +  v! .  (2.1.10) 

After  the  substitution  is  made,  the  resulting  equations  are  simplified.  The  Reynolds 
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stresses  are  the  terms  that  arise  in  the  Navier-Stokes  equations  when  Reynolds- 
averaging  the  momentum  equations.  The  term  is  given  by  T^e  =  pu'pi'y  The  Reynolds 
stresses  add  six  new  unknowns  into  the  equations  that  represent  the  transfer  of  mo¬ 
mentum  due  to  turbulent  fluctuation. 

Within  the  scope  of  RANS  turbulence  models  there  are  several  different  methods. 
The  least  complex  model  is  known  as  a  0-equation  model.  It  is  an  algebraic  model 
that  calculates  p  based  solely  on  local  flow  conditions.  A  1/2-equation  model  solves  an 
ODE  only  in  space  to  model  the  transport  of  turbulence.  A  1-equation  model  solves 
one  partial  differential  equation  (PDE)  of  the  transport  of  turbulence.  The  Spalart- 
Allmaras  (SA)  is  an  example  of  a  1-equation  model.  Two-equation  models  use  two 
turbulence  quantities  which  are  modeled  by  two  conservation  equations.  Examples  of 
2-equation  models  are  k  —  e  and  k  —  u. 

The  Boussinesq  gradient  transport  hypothesis,  also  known  as  the  eddy-viscosity 
hypothesis,  is  an  example  of  a  0-equation  model.  The  hypothesis  assumes  that  the 
fluxes  due  to  the  turbulent  fluctuations  could  be  solved  in  a  similar  manner  to  molecu¬ 
lar  transport  [4],  It  assumes  that  turbulent  shear  stress  depends  linearly  on  the  mean 
rate  of  strain.  These  transport  laws  can  then  be  defined  by  analogy  from  kinetic 
principles.  By  applying  this  hypothesis  to  the  Reynolds-averaged  stresses,  the  vis¬ 
cosity  term  is  replaced  by  a  sum  of  the  laminar  viscosity  and  the  turbulent  viscosity. 
Eddy-viscosity  is  not  valid  in  flows  with  sudden  change  of  mean  strain  rate,  signifi¬ 
cant  streamline  curvature,  rotation  and  stratification,  ducts  and  turbomachinary  and 
flows  with  boundary  layer  separation  and  reattachment. 

The  SA  model  solves  the  Reynolds-averaged  Navier-Stokes  equations  and  a  trans¬ 
port  equation  for  the  turbulence  model.  The  SA  model  has  the  benefit  of  being  local, 
numerically  forgiving,  and  rapid  convergence  to  steady  state  [29].  Reynolds-averaged 
equations  are  the  most  simplified  case  for  solving  turbulence  models,  if  not  the  most 
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widely  used.  The  comparatively  low  computational  cost  make  RANS  methods  appli¬ 
cable  to  engineering  applications;  however,  RANS  models  have  difficulty  predicting 
transition  to  turbulence.  [9] 

A  discussion  on  turbulence  modeling  would  be  incomplete  without  addressing 
DDES  [28].  Detached  eddy-simulation  (DES)  and  its  latest  form,  DDES  (Delayed 
Detached  eddy-simulation)  are  a  hybrid  of  RANS  and  LES.  DES  and  DDES  depend 
upon  a  defined  distance  to  the  wall.  Near  the  wall  the  models  act  as  a  RANS  model. 
When  the  distance  is  larger  than  a  predefined  distance,  the  model  acts  as  an  LES 
model.  DDES  has  modifications  to  ensure  that  RANS  turbulence  modeling  methods 
are  utilized  in  the  boundary  even  in  situations  with  thick  boundary  layers  as  well  as 
areas  of  shallow  separation.  The  need  for  DDES  is  motivated  by  the  grids  in  Figure 
3,  where  the  bottom  grids  would  initiate  a  switch  to  LES  methods  before  the  edge  of 
the  boundary  layer. 

2.1.3  Expected  Flow  Features 

It  is  important  to  determine  what  type  of  flow  will  result  in  the  simulation.  If  the 
expected  flow  regime  is  nearly  inviscid,  it  is  prudent  and  responsible  to  run  an  Euler 
solver  rather  than  a  turbulent  Navier-Stokes  solver  in  order  to  save  computation  time, 
effort  and  resources.  Flight  conditions  play  an  important  role,  as  around  Mach  0.8 
the  aircraft  or  airfoil  will  enter  into  the  transonic  regime.  Likewise,  at  very  high  Mach 
numbers  there  is  also  the  possibility  for  additional  thermophysical  considerations. 

The  angle  of  attack  will  play  a  large  factor  in  the  expected  results.  In  high  angle 
of  attack  (AoA)  flows,  highly  separated  turbulent  flow  is  to  be  expected.  In  the  case 
of  this  research,  the  SID  process  has  risen  out  of  the  need  to  model  the  highly  non¬ 
linear,  separated  flow  regime.  Therefore,  highly  separated  flow  is  expected  in  this 
research. 
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Figure  3.  Grids  in  a  boundary  layer.  Top  natural  DES;  left  ambiguous  spacing;  right 
LES.  Dotted  lines  mean  velocity.  S  is  the  boundary- layer  thickness  [27] 
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Vorticity  is  a  measure  of  the  rotation  that  exists  within  a  flow.  One  source  for 
this  rotation  is  viscosity.  In  the  Navier-Stokes  equations,  viscous  effects  arise  from 
the  viscous  stresses,  given  by  Equation  (2.1.4).  It  is  viscosity  that  creates  the  flow 
phenomena  necessary  for  flight,  such  as  boundary  layers.  Viscosity  also  introduces 
simulation  challenges  such  as  turbulence  and  separation.  For  a  high  fidelity  simula¬ 
tion,  solving  for  viscous  effects  are  very  important. 

Separation  occurs  as  the  boundary  layer  which  forms  on  an  airfoil  can  no  longer 
remain  attached  and  then  separates  from  the  airfoil.  The  point  at  which  separation 
occurs  is  dependent  on  several  factors,  including  angle  of  attack,  Mach  number  and 
Reynolds  number. 

2.2  Stability  and  Control 

Stability  and  control  are  fundamental  to  the  successful  flight  of  all  aircraft.  As  this 
project  aims  to  determine  these  important  parameters,  a  brief  introduction  regarding 
the  factors  to  be  used  will  be  presented  here. 

To  say  that  an  aircraft  is  stable  in  part  means  to  say  that  the  aircraft  is  statically 
stable.  It  is  important  that  when  an  aircraft  is  perturbed,  it  returns  to  a  trim  position. 
This  is  shown  in  Figure  4.  When  the  aircraft  is  trimmed,  the  static  stability  is 
measured  in  what  the  initial  reaction  of  the  system  is  to  a  perturbation.  If  there  is 
a  change  in  a,  such  as  a  gust,  and  there  is  no  aerodynamic  moment,  it  is  said  to  be 
neutrally  stable.  If  the  aircraft  returns  to  trim,  it  is  said  to  be  stable.  If  the  aircraft 
diverges  from  the  trim  conditions  it  is  said  to  be  unstable. 

While  static  stability  is  the  initial  tendency  of  a  system  to  respond  to  a  pertur¬ 
bation,  dynamic  stability  is  how  the  system  responds  to  the  perturbation  over  time. 
Dynamic  stability  is  shown  in  Figure  5.  Case  A  is  positive  static  stability  as  it  initially 
tends  towards  the  pre-displacement  state.  Case  A  is  also  dynamically  stable  as  over 
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Figure  4.  Stability  conditions  [31] 


time  it  returns  to  the  equilibrium  state.  Case  B  has  neutral  dynamic  stability  as  over 
time  it  does  not  tend  towards  equilibrium  nor  does  it  diverge.  Case  C  has  negative 
dynamic  stability  because  it  diverges  from  the  equilibrium  state. 

Insight  into  the  static  stability  of  an  aircraft  can  be  gained  through  the  stability 
coefficients  that  are  common  to  aeronautical  engineering  practices.  The  most  com¬ 
monly  sought  are  the  coefficient  of  lift  (Cl),  coefficient  of  drag  (Co),  and  the  pitch 
moment  coefficient  (Cm)-  These  coefficients  are  given  by  equation  (2.2.1) 
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2.3  System  Identification 


(2.2.1) 


System  Identification  is  one  of  three  main  problems  dealt  with  in  aircraft  dynamics 
and  control.  Referencing  Figure  6,  the  three  problems  are  [20]: 
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Figure  5.  Dynamic  stability  [3] 


1)  Simulation:  given  input  u  and  the  system  S,  find  the  output  y 

2)  Control:  given  the  system  S  and  the  output  y,  find  the  input  u 

3)  System  Identification:  given  the  input  u  and  output  y,  find  the  system  S 

The  methods  used  to  determine  the  system  are  referred  to  as  regression  methods. 
Regression  methods  are  a  statistical  technique  for  modeling  the  relationships  between 
variables.  For  example,  for  this  research  a  potential  model  could  be 

Cl  —  Ci  +  c^ol  +  c^cxQ  +  c^oi^Q  (2.3.1) 

where  cn  refers  to  a  constant  coefficient. 

Multivariate  orthogonal  functions  are  used  to  decorrelate  the  modeling  functions. 
Multivariate  orthogonal  functions  are  preferred  over  ordinary  least-squares  linear 
regression  techniques  as  it  is  difficult  to  identify  the  difference  between  correlated 
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system,  S 

Figure  6.  Aircraft  dynamics  and  control  [20] 


polynomials  obtained  through  ordinary  least-squares.  Once  the  regressors  have  been 
orthogonalized,  the  model  output  can  be  defined  as  [20] 

~  =  Pa  +  P  (2.3.2) 

where  z  is  a  vector  of  model  output,  a  is  a  vector  of  unknown  parameters,  P  is  a 
matrix  of  orthogonal  regressors  and  V  is  a  vector  of  measurement  errors. 

Put  simply,  SID  is  a  process  by  which  using  known  inputs  and  measured  outputs, 
a  complex  system  is  reduced  to  simplified  model.  In  this  work,  the  known  inputs 
are  parameters  such  as  a  and  Q.  The  measured  outputs  (5*)  will  be  the  stability 
coefficients,  Cl,  Cd,  and  Cm-  System  identification  is  a  method  by  which  to  cre¬ 
ate  a  mathematical  model  of  a  system  using  input  data.  This  method  with  CFD 
incorporated  is  outlined  in  Figure  7. 

There  are  three  main  steps  of  system  identification:  model  training,  model  vali¬ 
dation  and  model  prediction.  Model  training  is  the  act  of  creating  the  mathematical 
model  in  which  a  CFD  simulation  known  as  a  training  maneuver  is  conducted.  A 
training  maneuver  is  a  dynamic  computational  run  that  spans  a  desired  portion  of  the 
regressor  space.  The  dynamic  motion  of  the  training  maneuver  is  tailored  such  that 
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Figure  7.  The  system  identification  process  [18] 

the  aircraft  will  experience  a  useful  range  of  regressors  (a,  a,  Q,  Q,  etc)  values  that 
the  system  is  likely  to  experience  in  the  comparison  or  virtual  flight  test  maneuver 
stage. 

After  the  dynamic  computational  run  is  conducted,  the  reduced  order  mathemat¬ 
ical  model  is  fashioned  to  the  resulting  data.  This  mathematical  model  is  created  by 
examining  the  stability  characteristics  calculated  during  the  training  maneuver.  The 
model  will  predict  the  stability  characteristics  that  are  observed  in  a  real-world  flight 
test  maneuver.  There  are  several  different  methods  that  can  be  used  to  create  the 
reduced  order  mathematical  model.  Two  of  these  are  multivariate  polynomial  and  ra¬ 
dial  basis  functions.  Multivariate  polynomial  (MVP)  models  take  the  regressors  and 
raise  them  to  various  powers  along  with  coefficients  for  each  regressor  which  must  be 
uniquely  determined  for  a  given  system.  Radial  basis  functions  (RBFs)  are  based  on 
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the  radius  from  an  origin  or  center  and  generally  take  the  form  [25] 


V(x)  =  '}2wi(P(\\x  ~xi\\)  (2.3.3) 

i= 1 

where  y(x)  is  the  approximate  function,  wy  is  a  weighted  coefficient  and  Xi  is  a  moving 
center. 

Once  the  model  is  created,  the  SID  process  moves  to  Step  2,  model  validation. 
Model  validation  involves  demonstrating  that  the  reduced  mathematical  model  gives 
the  correct  coefficient  values.  The  reduced  model  is  applied  to  a  validation  maneuver 
and  compared  against  known  data  or  comparison  data  such  as  validated  CFD  results. 

If  the  model  passes  Step  2,  it  moves  on  to  Step  3,  model  prediction.  The  model  is 
then  used  to  predict  what  will  happen  to  an  aircraft  in  a  variety  of  maneuvers.  This 
point  in  the  process  is  also  known  as  “virtual  flight  testing”,  as  the  model  can  be 
used  to  “fly”  an  aircraft.  Using  the  developed  model,  time  and  money  will  be  saved 
by  reducing  the  amount  of  actual  flight  testing.  Virtual  flight  testing  can  indicate  to 
engineers  which  test  points  could  be  trouble  for  an  aircraft  and  which  ones  will  not 
be  an  issue. 

In  the  past,  wind  tunnel  testing  and  flight  testing  have  been  the  sole  inputs  to 
system  identification  methods.  A  variety  of  tools  have  been  developed  based  of  this 
type  of  input,  such  as  NASA’s  System  Identification  Programs  for  Aircraft  (SIDPAC) 
[23].  These  tools  focus  on  the  creation  of  the  reduced  model  and  generation  of  the 
resulting  data,  whereas  the  input  data  is  up  to  the  user.  It  is  only  recently  that  CFD 
has  been  explored  as  an  approach  for  system  identification. 
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2.4  Previous  Research 


This  section  will  look  into  the  research  that  has  previously  been  conducted  in  this 
field  and  will  consider  the  potential  shortcomings  of  that  research. 

2.4.1  Computational  Methods  for  Stability  and  Control  and  Early  At¬ 
tempts 

The  various  problems  discussed  in  Section  1.1  were  taken  from  a  NASA-sponsored 
symposium  of  S&C  and  CFD  engineers  known  as  Computational  Methods  for  Stabil¬ 
ity  and  Control  (COMSAC).  The  purpose  of  COMSAC  was  to  foster  the  discussion 
between  CFD  and  S&C  engineers,  fields  that  have  traditionally  has  been  considered 
two  separate  domains.  The  conclusions  drawn  from  the  symposium  include  [14]: 

•  The  inaccurate  prediction  of  aerodynamic  S&C  parameters  have  negative  ef¬ 
fects  on  the  life-cycle  costs  of  all  types  of  aircraft.  Generally,  the  inaccurate 
predictions  lead  developers  to  resort  to  “fly-and-try”  fixes. 

•  Improved  prediction  of  the  impact  of  separated  flow  on  the  aircraft  as  well  as 
the  character  of  the  separated  flow  should  be  a  priority. 

•  There  is  an  attitude  of  skepticism  in  using  CFD  for  aircraft  S&C  issues  in  both 
the  S&C  and  CFD  communities. 

•  The  success  of  advanced  CFD  methods  will  depend  on  the  demonstration  of 
success  for  generic  as  well  as  specific  aircraft  configurations. 

•  CFD  shortfalls  need  to  be  high  priority  targets  in  encouraging  the  use  of  CFD 
for  aircraft  S&C  tasks. 

•  There  was  immense  value  in  sharing  the  experiences  of  CFD  specialists  and 
S&C  specialists  at  the  symposium. 
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•  NASA  Ames  suggested  that  a  “brute  force”  method  be  utilized  [15]  to  fill  a  S&C 
database.  However,  it  was  determined  to  require  “30  different  angles-of-attack, 
20  different  Mach  numbers,  and  5  different  side-slip  angles,  each  for  a  number  of 
different  geometry  configurations  or  control  surface  deflections  [7]”.  Due  to  the 
large  computational  cost,  this  technique  would  need  to  rely  on  Euler  solutions 
which  was  deemed  undesirable. 

2.4.2  Training  Maneuvers 

The  training  maneuver  is  the  maneuver  which  will  provide  the  model  with  the 
known  inputs  and  the  measured  system  outputs.  In  this  case,  the  system  is  the  math¬ 
ematical  model  of  the  physical  system  utilized  in  the  CFD  solver.  It  is  very  important 
that  the  training  maneuver  be  able  to  accurately  excite  the  desired  phenomena  and 
be  general  enough  to  be  accurate  over  a  wide  range  of  scenarios.  Previous  research 
has  been  conducted  at  Air  Force  Institute  of  Technology  (AFIT)  [5],  United  States 
Air  Force  Academy  (USAFA)  [19]  [16],  the  Air  Force  Seek  Eagle  Office  (AFSEO)  [21] 
[12]  [13],  and  by  collaboration  [15]  [17]. 

McDaniel,  et  al.  [15],  examined  a  wide  variety  of  training  maneuvers  and  discussed 
how  well  the  resulting  SID  models  compared.  The  Erst  maneuver  examined  was  a 
pulse  plunge.  A  pulse  plunge  is  a  maneuver  designed  to  isolate  the  effects  of  a 
and  a  by  having  the  aircraft  vertically  drop  in  the  flow.  This  maneuver  highlights 
one  of  the  benefits  of  using  CFD  as  a  training  maneuver,  the  ability  to  simulate 
maneuvers  that  are  un-flyable,  meaning  pilots  in  actual  aircraft  cannot  command 
these  maneuvers.  Using  an  MVP  model,  the  results  were  reasonably  accurate  to  the 
validation  maneuver.  The  next  maneuver  is  the  Schroeder  sweep  plunge  maneuver, 
which  is  based  off  of  the  excitation  of  a  specified  bandwidth  through  the  summation 
of  discrete  frequencies.  Validated  against  a  DC  chirp  maneuver,  the  Schroeder  sweep 
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plunge  ROM  compared  accurately.  The  model  was  also  validated  against  a  pitch 
plunge  maneuver  where  the  model  less  accurately  followed  the  validation  maneuver! 
however,  it  was  noted  to  model  unsteady  behavior  well.  Then,  the  authors  trained  a 
model  using  the  DC  chirp  maneuver  previously  run  and  compared  it  against  sinusoidal 
pitch  motions.  The  resulting  models  were  comparable  for  a  maneuver  on  the  frequency 
of  2.0  Hz,  but  did  not  compare  as  well  for  1.0  Hz  or  for  0.5  Hz.  The  conclusion  was 
that  using  the  DC  chirp  input  the  MVP  model  could  model  the  system  response. 
However,  the  predictive  capability  was  limited  by  the  range  of  phenomena  excited  by 
the  training  maneuver. 

Dean,  et  al.  [13],  expanded  the  work  done  in  McDaniel,  et  al.  [15],  above  by 
then  comparing  the  results  against  solutions  obtained  by  Lockheed  Martin’s  6DOF 
simulator,  ATLAS.  Using  the  DC  chirp  maneuver  as  the  input,  the  SID  model  was 
compared  to  a  sinusoidal  pitching  maneuver  at  varying  frequencies.  At  1  Hz,  the 
model  had  trouble  predicting  CL  but  was  accurate  at  predicting  Cm-  As  the  frequency 
was  increased  to  2  Hz  and  3  Hz  the  model  was  able  to  better  predict  C-l-  However, 
it  was  not  able  to  correctly  predict  the  shape  of  the  hysteresis  curve.  It  should  be 
noted  that  the  authors  did  state  that  the  SID  models  take  seconds  to  compute  after 
the  CFD  training  maneuver  is  complete  and  are  able  to  offer  a  general  idea  of  the 
flight  characteristics,  whereas  ATLAS  data  is  only  available  after  significant  testing. 

Dean,  et  al.  [12],  again  expanded  upon  previous  work  by  comparing  SID  models 
created  from  CFD  against  a  flight  test  maneuver  conducted  by  a  F-16C  aircraft.  The 
maneuver  was  a  2.5-g  wind-up  turn  to  the  right.  The  authors  used  a  pitch  chirp 
maneuver  to  train  the  SID  model.  Lift  was  accurately  predicted,  while  drag  was 
observed  to  be  inaccurate  at  low  angles  of  attack.  The  cause  was  determined  to  be 
a  lack  inadequate  training  data  with  low  Q.  A  combined  training  maneuver  with 
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rotational  and  translational  motion  to  cover  the  whole  regressor  space  was  devised. 
The  results  for  drag  were  much  better  after  using  the  combined  training  maneuver. 

Some  work  has  also  been  done  by  Lillian,  et  al.  [21],  to  use  proper  orthogonal 
decomposition  to  recreate  the  pressure  distribution  on  a  wing  which  could  then  extract 
the  stability  coefficients. 

Jeans,  et  al.  [17],  examined  the  ability  of  DDES  simulations  to  predict  a  known 
non-linear  aerodynamic  effect  on  a  generic  fighter  geometry.  Then  the  group  used 
SID  techniques  to  created  ROMs  and  noted  how  well  the  dynamic  maneuvers  and 
the  ROMs  compared.  The  DDES  simulations  could  reliably  predict  the  instability 
within  a  given  frequency  range  and  the  SID  models  were  able  to  reproduce  constant 
frequency  maneuvers.  The  SID  techniques  were  also  found  to  have  worked  better  with 
higher  frequency  maneuvers.  The  authors  then  examined  how  the  SID  models  faired 
using  the  chirp  maneuver  as  the  training  maneuver.  The  group  found  that,  using 
the  SID  model  created  from  the  chirp  motion,  it  would  be  reasonable  to  extract  the 
rolling  moment  for  a  dynamic  roll  maneuver  with  frequencies  between  1.43  and  5.86 
Hz.  It  was  also  found  that  the  SID  models  indicate  unwanted  aerodynamic  behavior 
at  low  frequencies  which  would  alert  an  S&C  engineer  to  investigate  further. 

Dean,  et  al.  [11],  used  a  multitude  of  chirp  maneuvers  on  the  F-22  and  compared 
them  to  the  Lockheed  Martin  AVTEST  data.  The  training  maneuver  utilizes  a  chirp 
in  three  directions.  The  created  models  were  then  compared  against  a  variety  of 
maneuvers  such  as  static  data,  wind-up  turn  to  the  right,  and  a  pitch-up  stall  ma¬ 
neuver.  On  the  whole,  the  reduced  order  models  were  shown  to  compare  well  with 
the  validation  data. 

Jirasek,  et  al.  [19],  looked  again  to  study  which  types  of  training  maneuvers  pro¬ 
duced  desirable  results — specifically,  the  chirp,  DC  chirp,  spiral,  DC  spiral,  piecewise 
linear  spirals,  DC  piecewise  linear  spirals,  Schroeder,  plunging,  Fresnel  integral,  inte- 
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grated  DC  chirp,  integrated  spiral,  integrated  DC  spiral,  integrated  piecewise  linear 
spiral,  integrated  DC  piecewise  linear  spiral,  and  integrated  Schroeder  maneuvers. 
The  authors  found  that  the  spiral-  and  chirp-based  maneuvers  saw  better  lift  and 
pitching  moment  predictions.  Schroeder-based  models  resulted  in  favorable  predic¬ 
tions  at  lower  frequencies  and  unfavorable  predictions  at  higher  frequencies.  For  static 
data,  the  spiral  maneuver  gave  the  best  results.  Spiral-based  models  were  unable  to 
accurately  predict  sinusoidal  motion  while  the  chirp  and  Schroeder  maneuvers  were 
better  suited  to  that  task.  The  authors  state  that  aside  from  static  and  very  low 
frequency  data,  the  chirp  maneuvers  resulted  in  the  most  robust  and  reliable  models. 

2.4.3  Model  Generation 

There  has  also  been  research  into  which  type  of  models  to  use,  as  well  as  which 
regressors  are  important.  Two  models  considered  are  MVP  and  RBF. 

In  addition  to  the  training  maneuvers  discussed  above,  McDaniel,  et  al.  [15], 
looked  at  several  type  of  models,  specifically  multivariate  polynomial,  linear  regres¬ 
sion,  and  stepwise  regression  models.  The  results  showed  that  a  multivariate  polyno¬ 
mial  equation  better  compared  with  the  CFD  results  using  the  pitch  plunge  maneuver. 

Researchers  at  the  United  States  Air  Force  Academy  compared  MVP  models 
created  from  SIDPAC  software  and  radial  basis  function  models  to  wind  tunnel  data 
collected  on  site  on  the  X-31  [18].  The  conclusion  of  the  authors  was  that  the  SIDPAC 
and  RBF  models  were  equally  as  good  at  predicting  the  static  and  dynamic  pitching 
coefficients. 

Dean,  et  al.  [11],  examined  the  use  of  Kestrel  as  a  computational  resource  com¬ 
pared  to  Cobalt.  The  study  also  compared  MVP  and  RBF  models  created  using 
Kestrel  as  a  input  compared  against  Lockheed  Martin  performance  data.  Cl  and 
Cd  were  found  to  compare  well  for  both  models,  while  Cm  was  found  to  be  highly 
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inaccurate.  Minor  differences  between  the  models  can  be  seen  in  the  results,  however, 
the  differences  did  not  lead  the  authors  to  make  a  claim  as  to  which  model  is  more 
accurate. 

2.4.4  Reduced  Order  Modeling  of  Fighter  Aircraft 

Butler  [5]  began  to  look  into  the  objective  of  optimizing  training  maneuvers  based 
upon  a  series  of  what  he  called  ’’grid  metrics.”  Another  objective  was  to  improve  upon 
the  chirp  maneuver  such  that  it  was  better  suited  for  zero  Q  situations.  The  study 
was  conducted  using  a  full-scale  F-16  grid.  Two  problems  were  encountered  which  did 
not  enable  him  to  reach  any  conclusions  regarding  the  training  maneuver  parameters. 
First,  the  time  required  to  run  a  full-scale  F-16  grid  in  a  dynamic  maneuver  is  large, 
and  requires  a  vast  amount  of  computational  resources.  Running  several  maneuvers 
only  compounds  this  problem.  Thus,  the  work  of  Butler  was  only  brought  to  partial 
completion.  Second,  the  comparison  and  training  maneuvers  did  not  test  the  limits 
of  the  aircraft.  The  non-linear  coefficients  were  not  extremely  complicated,  making 
it  hard  for  there  to  be  much  of  a  difference  between  the  different  models.  Also,  it 
was  speculated  in  this  study  that  the  comparison  maneuvers  were  too  similar  to  the 
training  maneuvers  and  that  his  training  maneuvers  were  too  long,  giving  the  model 
more  than  ample  amount  of  data. 

Butler’s  work  was  invaluable  to  the  author  in  the  pursuit  of  this  research  and 
several  scripts  created  as  a  result  of  that  work  were  either  used  or  adopted  for  the 
author’s  use.  Also,  the  initial  RSPs  used  in  the  present  study  were  initially  proposed 
as  ‘grid  metrics’  in  Butler’s  work. 
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2.5  Kestrel 


The  CFD  software  chosen  for  use  in  this  research  is  the  product  of  the  Compu¬ 
tational  Research  and  Engineering  Acquisition  Tools  and  Environments  (CREATE) 
program  known  as  Kestrel.  Kestrel  began  as  an  initiative  from  the  Department  of 
Defense  High  Performance  Computing  Modernization  Program  (DoD  HPCMP)  [24], 
The  air  vehicles  portion  of  the  CREATE  program  is  known  as  CREATE-AV.  The 
job  of  the  CREATE-AV  team  was  to  determine  where  HPC  can  positively  affect  the 
acquisition  process.  It  was  determined  that  Kestrel,  a  virtual  fixed  wing  aircraft 
simulation  tool,  could  accomplish  this  goal. 

Kestrel  is  a  modular  program,  with  the  three  most  noteworthy  modules  or  com¬ 
ponents  being  the  Kestrel  Infrastructure  Executive  (KIE),  the  Kestrel  User  Interface 
(KUI),  and  the  flow  solver,  kAVUS,  a  derivative  of  the  Air  Vehicles  Unstructured 
Solver  (AVUS). 

KIE  is  the  executive  component  of  Kestrel.  It  is  the  job  of  KIE  to  push  data 
from  one  component  to  the  next.  KIE  however  does  not  ever  actually  utilize  the 
data  calculated  from  the  other  components.  KIE  is  event  based,  and  with  each  event, 
pushes  the  correct  data  pointer  to  the  component.  The  tasks  of  KIE  are  outlined 
below  [8]. 

a)  Read  the  XML  input  file  and  parse  the  data  out  to  each  component 

b)  Initiate  all  of  the  components. 

c)  Connect  all  data  consumers  with  producers. 

d)  Handle  the  publication  and  subscription  of  all  events. 

e)  Handle  all  component  exceptions. 

f)  Shut  down  all  components. 


27 


The  simplified  architecture  for  Kestrel  is  shown  in  Figure  8. 


Single  Executable  Additional 

Of  Modules  Executables 


Figure  8.  Kestrel  architecture  [24] 


KUI  is  the  interface  utilized  by  the  user  for  several  functions.  As  of  Kestrel  v3.0, 
KUI  has  6  modes  including  job  input,  boundary  condition,  mesh  manipulation,  post 
processing  and  preferences.  The  most  utilized  is  the  job  input  mode.  Here  the  user 
can  clearly  see  all  the  different  ‘knobs’  that  are  used  to  set  up  a  CFD  simulation  and 
does  not  have  to  interpret  the  meanings  of  various  entries  from  a  cryptic  job  file. 

The  most  important  component  for  this  work  is  kAVUS,  which  is  fundamentally 
a  finite- volume,  cell-centered,  first-order  in  space  and  time,  mesh-aligned  exact  Rie- 
mann  solver  of  Godunov.  kAVUS  is  modified  to  achieve  second-order  accuracy  in 
time  and  space.  Second-order  accuracy  in  space  is  patterned  after  van  Leers  MUSCL 
scheme  where  the  flow  state  is  assumed  to  vary  linearly  within  each  cell  [8].  First- 
and  second-  order  temporal  accuracy  is  achieved  via  the  unconditionally  stable  point- 
implicit  scheme  as  implemented  by  Tornaro  et  al  [30].  Also,  second-order  accurate 
viscous  terms  are  added  such  that  kAVUS  is  a  Navier-Stokes  solver.  For  turbulence 


modeling,  kAVUS  provides  the  following  options:  SA,  SA  +  DDES,  SARC,  SARC  + 
DDES,  Mentor,  Mentor  +  SST  and  Mentor  +  SST  +  DDES  [8], 
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III.  Methodology 


3.1  Overview 


The  methodology  of  this  thesis  is  organized  into  distinct  stages.  First,  the  grid 
to  be  used  for  the  remainder  of  the  research  was  created  and  validated.  The  initial 
RSPs  were  utilized  from  previous  research  conducted  by  Jed  Butler  [5].  These  RSPs 
will  be  assessed  against  a  known  robust  training  maneuver  and  known  low  quality 
training  maneuvers  in  order  to  determine  which  RSPs  provide  insight  into  the  training 
maneuver  development  process.  The  analysis  will  be  conducted  according  to  the 
overarching  process  depicted  in  Figure  9. 
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Figure  9.  Methodology  flow  chart 


First,  the  training  maneuvers  will  be  created  and  measured  using  the  regressor 
space  parameters.  Then,  the  training  maneuver  will  be  simulated.  The  resulting  co- 
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efficients  will  be  used  to  generate  the  reduced  order  models.  The  process  is  repeated 
for  each  training  maneuver;  then,  the  comparison  maneuvers  are  generated.  The  mo¬ 
tion  of  the  comparison  maneuvers  is  fed  into  the  ROMs  to  predict  the  comparison 
maneuver  coefficients.  That  prediction  is  compared  against  the  actual  comparison 
maneuver  CFD  coefficients.  The  metrics  R 2  and  normalized  root  mean  squared  devi¬ 
ation  (NRMSD)  are  then  used  to  gage  the  accuracy  of  the  predictions.  The  statistical 
analysis  program  JMP®is  then  used  to  examine  the  relationship  between  the  metrics 
and  the  RSPs.  All  computation  runs  will  be  conducted  at  Mach  0.5  at  10,000  ft. 

3.2  Computational  Resources  and  Hours 

This  research  was  conducted  using  two  computational  resources.  The  first  is  the 
Nordic  cluster  on-site  at  AFIT.  Nordic  is  a  12-node  system  comprised  of  16  processors 
per  node  with  two  of  the  processors  having  32  processors.  The  processors  are  2.4 
GHz  Opteron  with  64  (or  128)  GB  of  RAM  per  node  [1].  Jobs  are  submitted  via  the 
Portable  Batch  System  (PBS). 

The  second  computational  resource,  Raptor,  is  a  cluster  available  through  the  U.S. 
Air  Force  Research  Laboratory  (AFRL)  DSRC.  Raptor  has  2732  nodes  with  32  cores 
per  node.  Each  core  is  an  AMD  Opteron  64-bit  running  at  2.5  GHz  with  64  GB  of 
memory  [2], 


Table  1.  Computational  hours 


Cluster 

Total  CPU-Hours 

Grid  Study 

Raptor 

402919 

0 

Nordic 

14500 

14500 

Table  1  shows  the  number  of  computational  hours  used  on  each  cluster. 
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3.3  Time  Step  and  Grid  Density  Study 


The  geometry  chosen  for  this  project  is  a  NACA  64A010  airfoil.  The  decision 
to  utilize  a  2-D  case  was  made  in  order  to  reduce  computational  time  and  to  allow 
for  greater  simplicity  in  the  grid  generation  by  removing  difficulties  associated  with 
complex  geometries.  A  C-Grid  topology  was  chosen  for  ease  of  generation  which 
allowed  for  a  simple  way  to  modify  the  grid  refinement.  The  downside  of  a  C-Grid 
topology  is  that  there  are  very  refined  cells  where  there  is  no  need  for  them.  However, 
due  to  the  already  low  computational  expense  of  a  2-D  airfoil  case,  the  additional 
computational  expense  of  these  ‘extra’  cells  is  not  a  significant  consideration.  The 
coarse  grid  can  be  seen  in  Figure  10.  The  approach  taken  to  conduct  the  time-step 
and  grid  density  study  is  discussed  in  Cummings  et  al  [10]. 
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several  time  steps  are  also  chosen  using  varied  level  of  temporal  refinement.  The  rule 
of  thumb  of  At*  =  0.01  is  a  convenient  medium  or  fine  level  of  refinement,  where 


At* 


Atf/oo 

L 


(3.3.1) 


Then,  the  multiple  grids  are  run  in  a  systematic  way  using  the  different  time  steps. 
From  the  resulting  data,  a  spectral  analysis  is  conducted.  Using  the  normal  force 
coefficient,  a  single-sided  amplitude  spectrum  is  calculated  from  the  frequencies  ob¬ 
served  in  the  simulation  results.  The  dominant  frequency  is  pulled  from  this  analysis, 
and  the  Strouhal  number  is  calculated  by 


(3.3.2) 


where  /  is  the  frequency,  L  is  a  characteristic  length  and  U, A  is  the  flow  velocity.  In 
this  study,  L  is  taken  as  the  chord  length  of  1  ft.  The  wave  number  is  assumed  to 
be  the  inverse  of  the  Strouhal  number  and  is  plotted  versus  the  time  step.  Ideally,  a 
plot  like  Figure  11  will  result. 

Looking  at  Figure  11,  it  can  be  seen  that  the  fine  grid  and  the  medium  grid  con¬ 
verge  upon  the  same  Strouhal  number  at  a  time  step  of  2.5e-05  seconds.  It  would  be 
reasonable  then  to  run  either  of  the  grids  at  this  time  step,  however  for  computational 
expense  considerations,  it  would  be  wise  to  run  with  the  lower  grid  density.  The  main 
takeaway  from  this  graph,  however,  is  the  combined  effect  that  grid  density  and  time 
step  had  on  the  solution.  If  time  step  were  only  varied  for  a  single  grid  density,  one 
may  find  a  converging  Strouhal  number,  but  in  the  case  of  the  coarse  grid  in  Figure 
11,  the  convergence  behavior  would  be  inconsistent  with  a  fine  grid.  If  one  were  to 
vary  grid  density  for  a  set  time  step,  one  may  observe  convergence,  again  for  an  in¬ 
correct  value  as  is  seen  by  the  second  largest  time  step  in  Figure  11.  The  coarse  and 
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Figure  11.  Time-step  and  grid  convergence  study  for  NACA  6512  airfoil  [10] 


Table  2.  Grids  used  for  NACA  64A010  airfoil 


Grid  Density 

No.  of  Points 

No.  of  Cells 

No.  of  Faces  on  Surface 

Coarse 

124,125 

123,504 

248 

Medium 

298,850 

297,804 

498 

Fine 

802,575 

801,054 

498 

medium  grid  converge  on  a  Stronhal  number,  and  a  researcher  may  decide  to  run  at 
that  larger  time  step  on  the  coarse  grid  which  we  can  clearly  see  is  insufficient,  as 
proven  by  the  data. 

Table  2  shows  summary  data  for  the  three  grids  created  for  the  grid  density  study. 
Table  3  shows  the  different  time  steps  and  corresponding  At*  and  required  number  of 
iterations  for  0.4  seconds  of  simulation  time,  solutions  were  obtained  on  each  of  the 
three  grid  densities  from  Table  3  to  the  required  number  of  iterations  for  0.4  seconds 
of  runtime.  Then,  a  spectral  analysis  was  calculated  from  the  resulting  data  in  order 
to  fold  the  dominant  frequency  in  the  flow. 

The  spectral  analysis  was  conducted  using  a  fast  Fourier  transform  on  the  fluctu- 
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Table  3.  Time  steps  used  for  NACA  64A010  airfoil 

Time  Step  (s)  At*  Number  of  Iterations  for  0.4  sec 


5.0e-5 

0.0260 

8,000 

2.5e-5 

0.0134 

16,000 

1.25e-5 

0.0060 

32,000 

6.25e-6 

0.0030 

64,000 

ation  of  the  coefficient  of  normal  force  from  the  time-averaged  coefficient  of  normal 
force.  The  MATLAB®code  used  for  the  analysis  has  been  placed  in  Appendix  A.l  for 
the  reader’s  reference 

The  results  of  the  time  step  and  grid  density  study  are  discussed  in  Section  4.1. 

3.4  Initial  Regressor  Space  Parameters 

The  basis  for  this  research  started  with  the  “metrics”  provided  in  Butler  [5]  which 
have  been  re-designated  as  regressor  space  parameters  (RSPs)  in  this  work.  The 
RSPs  are  presented  below  and  are  calculated  using  the  MATLAB®script  contained  in 
Appendix  A. 3. 

The  goal  of  the  RSPs  is  to  quantify  how  the  training  maneuver  covers  the  regressor 
space.  For  this  research  the  regressor  space  is  considered  to  be  the  Q  values  between 
±75  (deg/s)  and  AoA  values  between  -5  and  20  degrees. 

The  RSPs  require  that  the  regressor  space  be  discretized  for  the  calculations.  For 
some  of  the  RSPs,  the  discretization  refinement  could  have  a  large  effect  on  the  results. 
The  question  is  how  fine  to  discretize  the  regressor  space.  While  pinpointing  an  exact 
number  is  not  imperative  (at  the  appropriate  range  the  same  trends  should  be  seen 
in  the  data)  finding  the  appropriate  range  is  important.  A  large  discretization  of  Q 
would  result  in  inflated  low  or  zero  Q  metrics.  However,  too  small  a  discretization 
would  exclude  values  from  zero  Q  calculations  that  are  practically  zero.  A  value  of 
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Q  and  AoA  discretization  of  0.1  was  recommended  by  Jed  Butler  [5]  and  is  used  in 
the  present  work  as  well.  The  first  RSP  is  an  overall  look  at  how  well  the  maneuver 
covers  the  desired  regressor  space. 


RSP] 


ft  of  cells  with  data  points 
Total  A  of  cells 


(3.4.1) 


The  second  metric  looks  to  measure  the  extent  to  which  the  extrema  of  the  regressor 
space  are  captured  by  the  maneuver. 


RSP2 


ft  of  cells  on  boundary  with  data  points 
Total  ft  of  cells  on  boundary 


(3.4.2) 


The  first  two  R.SPs  examine  how  many  points  total  the  TM  covers.  However,  it 
could  be  important  to  know  how  well  the  data  is  spread  throughout  the  regressor 
space.  This  is  intended  to  be  captured  by  Equation  (3.4.3).  In  this  RSP,  each  column 
of  discretized  AoA  values  is  counted  and  then  normalized  by  the  column  with  the 
maximum  number  of  points  in  it.  Then,  the  standard  deviation  is  taken  across  all 
the  columns.  In  an  evenly  spaced  training  maneuver,  the  standard  deviation  would 
be  large,  as  all  columns  would  have  the  same  number  of  points.  The  column  and  row 
descriptors  refer  to  the  format  presented  later  on. 


(  #  cells  with  data  points  i  ^ 

Total  #  of  cells  1 AoA  columnl 

RSP3  =  std 

#  cells  with  data  points  i 

Total  #  of  cells  1 AoA  column2 

max  value 

#  cells  with  data  points  i 
\  Total  #  of  cells  ' last  AoA  column  ) 
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A  similar  measure  of  evenness  is  calculated  using  the  rows  of  Q  values  instead  of 


columns  of  AoA  values.  This  is  shown  in  Equation  (3.4.4). 


(  #  cells  with  data  points  |  \ 

Total  #  of  cells  \Qrowl 

RSP4  =  std 

#  cells  with  data  points  i 

Total  #  of  cells  '  Q  row ^ 

4-  max  value 

#  cells  with  data  points  \ 

\  Total  #  of  cells  1 last  Q  row  ) 

(3.4.4) 


The  next  three  RSPs  focus  on  the  region  where  there  is  zero  pitch  rate.  The  Erst 
metric  shown  in  Equation  (3.4.5)  simply  takes  the  percentage  of  AoA  cells  with  zero 
Q  rate. 


ft  of  cells  at  0  Q  with  data  points 
R.o  f  5  =  - - — — — : — ; — - — - — 77 -  „  -  x  100 


(3.4.5) 


Total  ft  of  cells  at  0  Q 
The  next  RSP  shown  in  Equation  (3.4.6)  attempts  to  give  more  weight  to  training 
maneuvers  having  more  data  points  in  zero  Q  states.  This  hopes  to  allow  the  solution  a 
chance  to  dampen  out  dynamic  effects  and  gives  a  more  accurate  steady  state  solution. 
However,  it  is  possible  to  inflate  this  number  by  providing  a  large  amount  of  zero  Q 
data  that  does  not  add  any  more  accuracy  to  the  model.  This  RSP  is  also  influenced 
by  the  timestep  of  the  simulation. 


RSP*  = 


T otal  A  of  data  points  at  0  Q 
Total  ft  of  cells  atOQ 


(3.4.6) 


Similar  to  the  standard  deviation  RSPs  above,  this  next  RSP,  shown  in  Equation 
(3.4.7),  takes  the  standard  deviation  but  only  for  AoA  values  with  zero  Q. 


RSP7  =  std 


celli,  cell2,  cells,  ...last  cell 
Max  Value 


(3.4.7) 


@0Q 


The  next  set  of  three  RSPs  looks  to  capture  the  amount  of  low  Q  data.  Low  Q 
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data  can  be  defined  in  many  ways.  In  previous  research,  Bulter  [5]  looked  at  how 
different  definitions  affected  the  output  of  these  RSPs.  From  his  research,  common 
results  were  seen  when  the  cutoff  was  5  deg/s,  10  deg/s  and  20  deg/s  whereas  50 
deg/s  seemed  to  show  slightly  different  results.  Based  on  that  trend,  the  cutoff  for 
this  research  will  be  20  deg/s.  The  RSP  for  unique  cells  with  low  Q  is  shown  in 
Equation  (3.4.8). 


ff  of  cells  at  low  (not  0)  Q  with  data  points 
KoPr  = - — - - — - — - — — - - - ; - — — — -  x  100 


(3.4.8) 


Total  #  of  cells  at  low  ( not  0)  Q 
The  RSP  that  takes  into  account  all  the  data  points  at  low  Q  is  shown  in  Equation 
(3.4.9) 


RSP*  = 


Total  #  of  data  points  at  low  ( not  0)  Q 


(3.4.9) 


Total  ff  of  cells  at  low  ( not  0)  Q 
The  last  of  the  low  Q  RSPs  looks  at  the  standard  deviation  in  a  similar  fashion 
to  the  equations  prior.  This  RSP  is  shown  in  Equation  (3.4.10). 


RSP\o  =  std 


celli ,  cell 2,  cell3,  ...last  cell\  (3  4  10) 

Max  Value  )@iow(noto)Q 

The  last  set  of  the  initial  RSPs  look  at  the  opposite  side  of  the  20  deg/s  cutoff, 
and  are  considered  ‘high’  Q.  These  RSPs  are  computed  using  the  same  process  as  the 
last  three  RSPs,  but  using  the  high  Q  instead  of  the  low  Q  values.  These  three  RSPs 
are  shown  in  Equations  (3.4.11)  -  (3.4.13). 


ff  of  cells  at  high  Q  with  data  points 
11  Total  #  of  cells  at  high  Q 


(3.4.11) 


RSP*  = 


T otal  #  of  data  points  at  high  Q 
Total  #  of  cells  at  high  Q 


(3.4.12) 
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RSPit  =  std 


celli,  cell-2,  cell3,  ...last  cell 
M ax  V alue 


(3.4.13) 


@  high  Q 


3.5  Training  Maneuvers 


Unlike  the  previous  research  [5],  the  aim  of  this  project  is  not  to  test  out  new  train¬ 
ing  maneuvers.  The  focus  of  this  work  is  the  RSPs;  therefore,  the  training  maneuvers 
created  for  this  work  are  simpler  and  fewer  in  number.  The  standard  practice  ‘best’ 
training  maneuver,  a  DC  chirp  [19],  was  selected  to  be  the  baseline  maneuver.  Two 
other  training  maneuvers  were  created  to  be  the  ‘bad’  training  maneuvers,  designed 
to  highlight  the  weak  points  in  the  initial  RSPs  and  also  the  strong  points. 

Maneuvers  that  are  sinusoidal  in  nature  can  be  created  using  the  same  formula 
as  the  Kestrel  software,  and  then  can  be  input  into  Kestrel  as  arbitrary  motion  hies. 
For  the  purposes  of  this  discussion,  the  terms  will  be  left  in  the  form  used  by  Kestrel. 
‘S’  is  given  as  either  the  roll,  pitch  or  yaw  angle,  depending  on  the  axis  of  prescribed 
sinusoidal  motion. 


S(t )  =  s(t)cos  [2n(f3ft1+Xf  +  fit  +  <h/360)]  (3.5.1) 


where 


s(i)  =  (3ai1+Xa 

_  02  —  ai 

Pa 


( tmax  t0)A“ 


t=(t-  t0) 

fl-fl 


Pf  = 


(tmax  tg)  f 


(3.5.2) 

(3.5.3) 

(3.5.4) 

(3.5.5) 


where  the  variable  t  is  actual  time,  tmax  is  the  maximum  time  of  the  motion,  to  is  the 
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initial  time  of  the  motion,  a2  and  /2  are  the  final  amplitude  and  frequency,  whereas 
a  1  and  f\  are  the  initial  values.  A /  and  Aa  are  the  shift  parameters  that  define  how 
the  motion  will  be  computed  from  the  initial  to  the  final  values.  A  value  of  one  in 
either  A  parameter  represents  a  linear  change.  <f>  is  the  phase  shift  given  in  degrees 
[8]. 

Using  Matlab®,  these  equations  can  be  calculated  locally,  and  then  the  RSPs 
discussed  in  Section  3.4  can  be  calculated.  From  there  the  maneuver  is  transferred 
into  an  arbitrary  motion  hie  as  described  by  the  Kestrel  User  Guide  [8].  The  user 
guide  specifies  that  there  must  be  eight  header  lines  followed  by  13  columns.  The 
first  column  is  time.  The  next  three  columns  are  the  rotated  basis  vector  nx,  followed 
by  the  next  three  columns  of  ny  and  three  columns  of  nz.  The  last  three  columns  are 
the  location  of  the  center  of  rotation. 

Three  sinusoidal  training  maneuvers  were  created  using  Equations  (3.5.1)  -  (3.5.5). 
The  inputs  for  these  training  maneuvers  can  be  seen  in  Table  4.  It  was  quickly 
noticed  that  using  an  arbitrary  motion  hie  rather  than  the  Kestrel  inputs  significantly 
increased  the  computational  time  per  iteration,  training  maneuver  2,  comparison 
maneuver  la,  comparison  maneuver  3  and  comparison  maneuver  4  were  computed 
using  an  arbitrary  motion  hie,  the  rest  of  the  maneuvers  were  run  using  the  inputs  in 
KUI.  For  the  sinusoidal  maneuvers  such  as  training  maneuver  1,  training  maneuver 
3  and  comparison  maneuver  2,  these  inputs  are  identical,  comparison  maneuver  lb 
and  comparison  maneuver  lc  required  use  of  the  constant  rate  pitch-and-hold  motion 
in  Kestrel.  Comparison  maneuver  6  used  a  series  of  constant  rate  pitch  and  hold 
maneuvers.  Comparison  maneuver  2  uses  multiple  sinusoidal  inputs 

The  hrst  maneuver  is  the  best  practice  maneuver,  it  is  known  as  a  DC  Chirp 
maneuver.  The  goal  is  to  cover  a  large  range  of  the  regressor  space.  While  the  DC 
chirp  produces  the  most  robust  models,  it  is  known  to  have  inaccuracies  in  the  low 
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Table  4.  Training  maneuver  input  parameters 


Input 

TM1 

TM2 

TM3 

Aa 

1 

1 

1 

A  / 

1.9 

1 

1 

CL\ 

14.7 

12.5 

10 

d2 

0 

12.5 

10 

h 

0.75 

1 

0.75 

f'2 

1.6 

1 

0.5 

<F (degrees) 

-90 

-90 

-90 

AoAiidegrees) 

6.5 

7.5 

7.5 

Tirne(s) 

4 

1 

3 

and  zero  Q  ranges.  The  first  training  maneuver  (TM)  is  shown  in  Figure  12.  By 
covering  a  greater  portion  of  the  regressor  space,  it  is  desired  that  the  model  will  be 
more  accurate  with  a  greater  understanding  of  the  varying  flow  physics. 

This  training  maneuver  aims  to  span  a  large  amount  of  the  regressor  space.  To 
accomplish  that  task,  the  frequency  is  varied  as  is  the  amplitude  of  the  wave.  The 
result  is  seen  in  Figure  13. 

Figure  14  shows  a  3-D  plot  of  the  regressor  space.  Cells  with  greater  spikes  are  cells 
that  have  data  points  more  often  in  that  particular  cell.  For  visualization  purposes, 
the  cells  for  the  3-D  plots  are  considered  to  be  1  (deg  and  deg/s)  as  smaller  cells  make 
the  3-D  plots  impossible  to  see.  As  can  be  seen  with  Figure  14,  while  the  number  of 
points  centers  around  the  initial  AoA,  as  would  be  expected,  the  maneuver  covers  a 
variety  of  points  and  does  not  center  on  one  particular  region. 

TM2  and  TM3  are  the  maneuvers  designed  to  produce  less  than  optimal  results. 
TM2  is  a  simple  single  oscillation  sinusoid,  as  can  be  seen  in  Figure  15. 

TM2  reaches  the  boundary  of  the  regressor  space  during  the  maneuver,  hitting 
the  maximum  and  minimum  Q  and  AoA.  However,  it  does  not  make  any  effort  to 
cover  the  large  range  of  the  regressor  space  as  can  be  seen  in  Figure  16. 
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Figure  12.  Training  maneuver  1:  AoA  and  Q  vs  time 


AoA  (deg) 

Figure  13.  Training  maneuver  1:  Q  vs  AoA 
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Pitch  Rate  (Q) 


Number  of  Data  Points 


3000 


2000 


1000 


-50 

Pitch  Rate  (Q)  _ioo 


-10 


AoA 


20 


Figure  14.  Training  maneuver  1:  discretized  regressor  space 


Figure  15.  Training  maneuver  2:  AoA  and  Q  vs  time 
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Pitch  Rate  (Q) 


Figure  16.  Training  maneuver  2:  Q  vs  AoA 


The  last  plot  shown  in  Figure  17  shows  the  amount  of  hits  in  the  regressor  space. 
As  would  be  expected,  there  are  many  hits  along  the  outskirts  of  the  regressor  space, 
however  there  are  space  the  middle  of  the  regressor  space  has  no  data  points.  The 
number  of  data  points  overall  is  less  than  seen  in  Figure  14  for  TM1.  The  lower  of 
points  results  from  the  shorter  time  period  for  TM2  (1  second)  vs  the  time  period  for 
TM1  (4  seconds). 

The  last  maneuver,  TM3,  attempts  to  cover  slightly  more  of  the  regressor  space 
while  also  covering  more  low  Q  data.  However,  it  is  also  an  aim  of  this  training 
maneuver  to  miss  the  boundaries  of  the  regressor  space  entirely  as  can  be  seen  in 
Figure  18. 

The  emphasis  around  the  low  Q  threshold  of  20  deg/s  can  be  seen  in  Figure  18 
and  Figure  19.  The  amount  of  hits  in  the  regressor  space  can  be  seen  in  Figure  20. 
While  TM3  shows  a  greater  spread  than  TM2,  it  still  does  not  cover  as  large  a  range 
as  TM1. 
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Number  of  Data  Points 


Figure  17.  Training  maneuver  2:  discretized  regressor  space 


Figure  18.  Training  maneuver  3:  AoA  and  Q  vs  time 
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Pitch  Rate  (Q) 


Number  of  Data  Points  Pitch  Rate  (Q) 


AoA  (deg) 

Figure  19.  Training  maneuver  3:  Q  vs  AoA 


Pitch  Rate  (Q)  -ioo  _io 


Figure  20.  Training  maneuver  3:  discretized  regressor  space 
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Both  TM2  and  TM3  are  distinguishable  from  TM1  in  that  the  3-D  plots  are  much 
more  sharply  divided  than  the  3-D  plot  for  TM1.  The  sharp  contrast  for  TM2  and 
TM3  are  because  of  the  single  sweep  of  the  data  regime.  TM2  and  TM3  make  one 
pass  and  there  is  not  much  overlap  in  the  points. 

3.6  Comparison  Maneuvers 

The  purpose  of  the  comparison  maneuvers  is  to  provide  the  reduced  order  mod¬ 
els  generated  from  the  training  maneuvers  a  set  of  comparison  data.  That  is,  the 
comparison  maneuvers  are  simulated  using  Kestrel.  Instead  of  creating  reduced  order 
models  from  these  maneuvers,  the  results  of  the  comparison  maneuver  simulations 
are  used  as  a  basis  of  comparison  for  the  ROMs. 

The  first  comparison  maneuver,  abbreviated  as  COM1,  is  simple.  COM1  is  a 
linear  sweep  across  the  AoA  regime.  The  goal  of  these  maneuvers  is  to  easily  see 
how  well  the  models  can  recreate  the  entire  a  range  at  a  set  Q  value.  Therefore, 
COM1  actually  consists  of  three  different  maneuvers:  COMla,  COMlb  and  COMlc. 
COMla  is  set  to  a  Q  value  of  25  deg/s,  COMlb  has  a  Q  value  of  10  deg/s  and  COMlc 
has  a  Q  value  of  50  deg/s.  a  and  Q  for  the  3  maneuvers  can  be  seen  in  Figures  21  - 

23. 

The  goal  of  COM1  is  to  see  if  the  ROMs  can  produce  accurate  O  vs  a  curve 
and  Cd  arid  Cm  curves.  These  are  important  curves  which  are  the  basis  of  much  of 
aircraft  design. 

The  next  maneuver  known  as  COM2  is  more  complicated.  This  maneuver  at¬ 
tempts  to  cover  a  very  large  range  of  the  regressor  space  and  can  be  seen  in  Figure 

24.  While  going  through  the  oscillations  shown  in  Figure  24,  the  maneuver  also  covers 
a  large  range  of  Q  values,  which  can  be  seen  in  Figures  25  and  26. 

The  next  two  maneuvers  strive  to  go  through  the  AoA  region  for  a  range  of  Q 
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Figure  21.  Comparison  maneuver  la  (25  deg/s):  AoA  and  Q  vs  ti 


Figure  22.  Comparison  maneuver  lb  (10  deg/s):  AoA  and  Q  vs  time 
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Pitch  Rate  (Q)  5  Pitch  Rate  (Q) 


Figure  23.  Comparison  maneuver  lc  (50  deg/s):  AoA  and  Q  vs  time 


Figure  24.  Comparison  maneuver  2:  AoA  vs  time 
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Pitch  Rate  (Q) 


Figure  25.  Comparison  maneuver  2:  AoA  and  Q  vs  time 


AoA  (deg) 

Figure  26.  Comparison  maneuver  2:  Q  vs  AoA 
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Pitch  Rate  (Q) 


values.  COM3  uses  an  exponential  change  in  angle  of  attack  in  order  to  provide  a 
range  of  AoA  and  Q  values.  COM3  can  be  seen  in  Figure  27  and  28. 


g 
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Figure  27.  Comparison  maneuver  3:  AoA  and  Q  vs  time 


While  COM3  covers  the  range  of  positive  Q-values,  COM4  covers  the  negative 
Q-values  and  can  be  seen  in  Figures  29  and  30. 

The  final  maneuver,  known  as  COM5,  is  a  series  of  constant  rate  pitch  and  hold 
maneuvers.  The  goal  for  COM3  is  to  show  the  weakness  of  the  maneuvers  at  0  Q 
values  at  several  angles  of  attack.  The  maneuver  can  be  seen  in  Figure  31. 


3.7  Creating  the  Reduced  Order  Models 

The  reduced  order  models  are  created  using  methods  discussed  in  2.3.  The  soft¬ 
ware  used  to  conduct  this  analysis  is  the  System  Identification  Programs  for  Aircraft 
(SIDPAC)  produced  by  NASA  Langley.  For  the  purposes  of  this  research,  it  was  de¬ 
cided  to  not  use  the  built-in  post-processing  system  identification  tools  implemented 
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Pitch  Rate  (Q) 


AoA  (deg) 

Figure  28.  Comparison  maneuver  3:  Q  vs  AoA 


Figure  29.  Comparison  maneuver  4:  AoA  and  Q  vs  time 
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Pitch  Rate  (Q) 


Pitch  Rate  (Q) 


AoA  (deg) 

Figure  30.  Comparison  maneuver  4:  Q  vs  AoA 


Figure  31.  Comparison  maneuver  5:  AoA  vs  time 
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by  Kestrel.  This  was  done  to  ensure  greater  control  and  understanding  of  model  cre¬ 
ation.  The  SIDPAC  software  is  a  collection  of  MATLAB®functions,  and  the  usefulness 
of  each  one  is  dependent  upon  the  need  of  the  user.  For  this  research  the  code  mof.m 
will  be  used.  This  file  will  map  a  multivariate  polynomial  model  to  the  input  data. 
The  file  does  require  some  input  from  the  user.  The  script  used  to  call  mof.m  and 
input  the  necessary  parameters  is  shown  in  Appendix  A. 2. 

The  first  input  is  the  variable  nord ,  which  is  described  as  a  vector  of  maximum 
independent  variable  orders.  Each  input  variable,  a  and  Q,  is  given  a  nord  value. 
This  value  refers  to  maximum  order  of  each  independent  variable  in  the  final  reduced 
order  model.  If  a  4  is  entered  for  a,  the  largest  exponent  on  a  is  4  in  the  ROM. 
However  another  term  maxord  restricts  the  total  number  of  exponents  for  a  given 
term.  If  maxord  is  defined  to  be  5,  there  could  be  a  term  with  a4Q  but  not  a4Q2  as 
4+2=6  which  is  greater  than  5.  The  final  option  is  maxopt  which  sets  the  maximum 
number  of  terms  in  the  final  polynomial  model. 

For  this  work,  the  maxord  is  selected  to  be  5,  as  additional  terms  were  seen  to 
increase  computational  time  significantly  while  not  improving  the  ROM.  The  pro¬ 
gram  orders  the  terms  by  significance  to  the  model  and  ends  where  there  is  limited 
additional  benefit  to  the  model  based  on  R2  and  predicted  square  error  (PSE).  From 
experience  in  this  research,  the  maximum  value  was  seen  to  be  around  14  in  most 
cases.  However,  the  maximum  number  of  terms  was  not  limited  by  user  input  in  this 
research.  An  example  of  the  mof.m  output  is  shown  in  Figure  32.  The  red  dot  in 
the  figure  represents  the  programs  final  choice  for  the  number  of  terms  in  the  model. 
This  example  is  the  output  fit  for  TM1  Cl  data.  The  R2  value  in  this  chart  is  the 
R2  for  the  model  fit  to  the  raw  CFD  data.  It  is  lower  than  the  R2  values  shown  in 
Chapter  IV  due  to  the  fluctuations  in  the  flow. 

For  each  training  maneuver,  each  regressor  will  have  a  nord  value  of  5,  and  maxord 
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Prediction  Error 


Figure  32.  Model  termination  determination  by  SIDPAC  software 
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will  be  defined  as  5.  Then  a  model  will  be  generated  using  the  mof.m  program.  The 
program  will  calculate  the  best  model  given  the  input  data  and  input  parameters. 
That  reduced  order  model  will  then  be  used  going  forward. 

In  the  highly  constrained  circumstances  of  CFD  flow,  a  is  assumed  to  be  equal 
to  Q,  as  any  other  motion  will  not  be  introduced.  However  Q  could  also  be  a  factor 
for  the  model.  In  order  to  provide  Q  as  a  regressor,  the  deriv.m  function  of  SIDPAC 
will  be  used.  This  function  calculates  a  locally  smoothed  numerical  time  derivative. 
However,  Q  was  not  selected  as  a  regressor  in  any  of  the  final  reduced  order  models. 

3.8  Analysis  of  Results 

Once  the  reduced  order  models  have  been  created,  they  will  be  compared  to  full 
CFD  predictions  of  the  comparison  maneuver  coefficients  Cl,  Co  and  Cm-  Because 
of  the  unsteady  nature  of  the  flow,  a  moving  average  is  computed  for  the  coefficients. 
The  metrics  for  each  ROM,  as  compared  to  the  comparison  maneuver  moving-average 
coefficients,  are  calculated  using  the  following  relations: 


z  =  meaniObs) 

(3.8.1) 

SSR  =  ( Model  -  z )2 

(3.8.2) 

SSE  =  ( Obs  -  Model)2 

(3.8.3) 

SST  =  SSR  +  SSE 

(3.8.4) 

9  SSR 

(3.8.5) 

“  SST 

where  R2  is  the  coefficient  of  determination,  z  is  the  mean  of  the  observed  (CFD 
coefficient)  data,  Model  refers  to  the  ROM  predicted  coefficient,  SSR  is  the  regression 
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sum  of  squares,  SSE  is  sum  of  squares  of  residuals  and  SST  is  total  sum  of  squares. 


NRMSD 


^2  (Obs— Model)2 
#  of  observations 


max  (Obs )  —  min(Obs ) 


(3.8.6) 


NRMSD  is  the  normalized  root  mean  squared  deviation.  Each  ROM,  as  derived 
from  the  various  TMs,  will  have  these  two  metrics  associated  with  them  for  each 
comparison  maneuver. 

Two  approaches  will  be  used  to  analyze  the  results.  The  first  is  a  more  visual 
approach.  An  examination  of  the  RSPs  and  the  output  metrics  for  each  maneuver 
will  be  compared  along  with  plots  of  the  predicted  and  simulated  coefficients  to  get  a 
sense  of  which  RSPs  show  a  trend  in  the  final  data.  The  second  approach,  to  add  more 
scientific  rigor,  will  incorporate  the  statistical  program  JMP®to  investigate  which 
inputs  (RSPs)  have  an  effect  on  the  outputs  (R2  and  NRMSD)  for  each  coefficient 
using  “contrasting”  methods.  It  is  important  to  note  the  combined  effectiveness 
gained  by  using  these  two  methods.  While  JMP®will  be  able  to  provide  insight 
into  the  generic  trend  of  the  data,  it  will  be  the  onus  of  the  author  to  correlate 
trends  within  the  data  subsets  such  as  low  Q  maneuvers  or  even  low  Q  sections  of 
maneuvers. 
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IV.  Results 


4.1  Time  Step  and  Grid  Density  Study 

Using  the  methodology  discussed  in  Section  3.3,  the  time  step  and  grid  density 
study  was  conducted,  Figure  33  shows  the  wave  number  converges  between  the  fine 
and  the  medium  grid  at  a  time  step  of  2.5  x  10~5  seconds.  However,  at  this  point  the 
two  grids  are  still  decreasing  in  wave  number  with  respect  to  decreasing  time  step. 
At  a  time  step  of  1.25  x  10~5  seconds,  both  grids  are  stationary  with  respect  to  wave 
number  and,  are  therefore,  considered  to  be  converged  with  respect  to  time.  At  this 
point,  it  can  be  said  that  the  solution  is  grid  and  time  step  independent.  At*  =  0.01 
for  these  flow  conditions  is  plotted  for  reference.  It  is  desirable  to  run  on  the  smallest 
grid  possible  to  reduce  computation  expense;  therefore,  the  medium  grid  at  a  time 
step  of  1.25  x  10~5  seconds  was  chosen.  Figure  34  is  a  graph  of  the  time-averaged 
value  for  normal  force  coefficient.  While  not  a  great  indicator  of  grid  convergence, 
it  does  serve  as  a  sanity  check.  The  medium  grid  and  the  fine  grid  nearly  converge 
upon  a  solution  at  the  1.25  x  10-5  second  time  step. 

However,  after  the  medium  grid  was  selected  and  the  training  maneuver  simula¬ 
tions  began,  the  grid  required  additional  refinement.  Additional  cells  were  needed 
approximately  0.1  ft  above  and  below  the  trailing  edge.  The  final  grid  can  be  seen  in 
Figure  35  and  the  final  grid  statistics  are  given  in  Table  5. 

Table  5.  Final  grid  used  for  NACA  64A010  airfoil 

Grid  Number  of  Points  Number  of  Cells  Number  of  Faces  on  Surface 
Final  629,378  628,008  498 
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Time  Average  normal  force  coefficient  Wave  Number  (1/St) 


Figure  33.  Grid  sensitivity  study  results:  wave  number 


Figure  34.  Grid  sensitivity  study  results:  time-averaged  normal  force 
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Figure  35.  Final  computational  grid 


4.2  Initial  RSPs 

The  initial  RSPs  were  computed  for  the  three  training  maneuvers  and  the  results 
are  shown  in  Table  6  and  Table  7.  The  first  thing  to  note  is  the  addition  of  a  TM4.  It 
was  observed  early  on  that  it  would  help  the  statistical  analysis  if  there  was  another 
point  of  comparison  for  the  RSPs.  Since  the  COM2  CFD  results  were  readily  available 
and  COM2  was  very  similar  in  behavior  to  a  training  maneuver,  the  results  of  COM2 
were  used  to  create  additional  models  and  labeled  as  TM4.  However,  TM4  is  not 
compared  to  COM2  as  that  would  confound  the  results. 

Looking  at  the  RSPs  can  be  useful  for  the  purpose  of  gaining  a  better  idea  of 
what  the  maneuver  is  doing  before  looking  at  the  output  metrics.  From  the  first 
seven  RSPs  it  can  be  seen  that  TM4  and  TM1  cover  a  larger  amount  of  the  regressor 
space,  while  TM2  and  TM3  do  not.  This  is  to  be  expected,  as  these  maneuvers  were 
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Table  6.  First  seven  RSPs 


TM 

Percent 

Whole  Boundary 

Std 

a 

on 

Q 

Percent 

0  Q  Data 
Per  cell 

Std 

TM1 

3.6691 

4.7483 

0.1951 

0.2444 

5.20 

0.9720 

0.1000 

TM2 

0.8851 

11.556 

0.1074 

0.1532 

0.80 

0.1280 

0.0892 

TM3 

0.6880 

0.0000 

0.1188 

0.1105 

0.80 

0.9889 

0.0671 

TM4 

5.0029 

15.103 

0.1605 

0.1862 

4.80 

0.8320 

0.1404 

Table  7. 

Last  six  RSPs  with  20  deg/s 

>  as  cutoff 

TM 

Below  20  deg/s 
Percent  Per  cell  std 

Above  20  deg/s 
Percent  Per  cell  std 

TM1 

4.9990 

0.9586 

0.1226 

3.1866 

0.8152 

0.1053 

TM2 

0.8120 

0.1311 

0.0891 

0.9115 

0.1874 

0.0903 

TM3 

1.3050 

1.3431 

0.0757 

0.4650 

0.3845 

0.0336 

TM4 

4.2240 

0.8766 

0.1472 

5.2843 

1.1357 

0.1404 

not  created  to  cover  a  wide  range.  With  these  high  values  it  would  be  expected  that 
TM1  and  TM4  would  correlate  to  an  overall  robust  model. 

The  next  RSP  is  the  percent  of  the  boundary.  TM2  and  TM4  both  have  large 
values  for  this  RSP.  This  result  may  correlate  to  favorable  results  for  coefficients  that 
are  more  linear  in  nature. 

The  third  and  fourth  RSPs  are  the  standard  deviation  in  AoA  and  Q.  A  high  stan¬ 
dard  deviation  would  imply  that  the  results  are  more  evenly  distributed  throughout 
the  regressor  space.  A  small  standard  deviation  would  imply  that  a  majority  of  the 
data  set  is  located  near  the  mean.  TM1  provides  the  largest  standard  deviation  in 
AoA  and  Q. 

The  next  three  RSPs  focus  on  the  zero  Q  data.  The  same  thoughts  above  for  the 
general  case  carry  down  into  this  data  subset.  TM1  hits  the  largest  percent  of  the 
dataset  as  well  as  the  largest  per  cell  count  and  one  of  the  largest  standard  deviations. 
Therefore,  it  is  expected  to  produce  the  best  results  in  a  zero  Q  scenario. 
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The  next  set  looks  at  the  low  Q  data  points.  Looking  at  these  RSPs,  TM1  again 
returns  the  best  values.  However,  for  the  high  Q  RSPs,  TM4  returns  the  best  values. 
Based  on  these  values  it  would  be  expected  that  TM1  models  will  return  the  best  low 
Q  data  while  TM4  models  will  return  the  best  high  Q  data. 

4.3  Model  Results 

This  section  discusses  the  results  gained  from  comparing  the  reduced  order  models 
to  the  results  gained  from  the  full  CFD  simulation  of  the  comparison  maneuvers. 
Again,  TM4  is  not  compared  against  COM2  as  that  would  just  be  a  measure  of  how 
well  the  model  fit  the  data  and  not  a  measure  of  the  predictive  capability  of  the 
model. 

4.3.1  Coefficient  of  Lift 

The  results  for  Cl  can  be  seen  in  Table  8.  Looking  at  the  table  there  are  some 
surprising  results.  From  the  RSPs  and  the  way  the  training  maneuvers  were  designed, 
TM2  and  TM3  were  not  expected  to  produce  very  good  results.  However  it  is  clear 
from  the  metrics,  and  the  results  that  this  is  not  the  case  for  coefficient  of  lift.  To 
refresh  the  memory  of  the  reader,  COMla  is  the  medium  Q  maneuver  at  25  deg/s, 
COMlb  is  the  low  Q  maneuver  at  10  deg/s  and  COMlc  is  the  high  Q  maneuver  at 
50  deg/s. 

First,  reference  Figure  36,  it  can  be  seen  that  the  models  all  follow  the  trend  of 
the  comparison  data  well.  At  first  glance,  it  would  appear  that  all  the  models  are  a 
good  prediction,  with  TM4  being  the  only  model  distinguishable  from  the  rest  of  the 
data.  Table  8  confirms  that  idea.  The  surprising  result  is  that  TM2  has  the  highest 
R2  and  the  lowest  NRMSD.  Following  that  the  best  model  is  TM3.  Linking  back 
to  the  RSPs  in  Tables  6  and  7,  TM2  is  middle  of  the  pack  on  everything  but  the 
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Table  8.  R 2  and  NRMSD  for  all  Cl  models  vs  a  moving  average 


COM 

TM1 

TM2 

TM3 

TM4 

R2 

NRMSD 

R2 

NRMSD 

R2 

NRMSD 

R2 

NRMSD 

COMla 

0.9920 

0.0254 

0.9941 

0.0211 

0.9936 

0.0216 

0.9891 

0.0271 

COMlb 

0.9926 

0.0237 

0.9939 

0.0211 

0.9918 

0.0256 

0.9860 

0.0305 

COMlc 

0.9919 

0.0227 

0.9905 

0.0245 

0.6674 

0.1068 

0.9917 

0.0219 

COM2 

0.9774 

0.0376 

0.9826 

0.0329 

0.8064 

0.1034 

- 

- 

COM3 

0.9950 

0.0197 

0.9941 

0.0212 

0.6487 

0.1707 

0.9947 

0.0407 

COM4 

0.9619 

0.0407 

0.9809 

0.0300 

0.7761 

0.0912 

0.9823 

0.0283 

COM5 

0.9901 

0.0245 

0.9904 

0.0244 

0.9852 

0.0318 

0.9822 

0.0318 

zero/low  Q  parameters  in  which  it  is  the  lowest,  therefore  there  is  no  clear  link  from 
the  metrics  to  the  RSPs  from  this  maneuver. 

Next,  reference  Figure  37.  The  results  for  this  case  are  very  similar  to  the  results 
for  COMla.  Looking  at  the  results  in  Figure  37  it  is  hard  to  distinguish  the  results 
and  all  the  models  agree  well.  Looking  at  the  metrics  in  Table  8  it  can  be  seen  that 
TM2  has  the  largest  R2  and  lowest  NRMSD.  Surprisingly,  despite  having  a  greater 
concentration  in  the  low  Q  area,  TM3  has  worse  values  than  TM1  if  only  by  a  small 
margin.  Because  of  this,  there  is  again  no  clear  link  between  the  output  metrics  and 
the  RSPs. 

The  last  constant-Q  a  sweep  results  are  shown  in  Figure  38.  This  high  Q  case  is 
where  we  start  to  see  some  separation  between  the  models.  For  this  maneuver  it  is 
easily  seen  that  TM3  is  wildly  off  point.  The  error  is  due  to  the  fact  that  TM3  does 
not  have  a  large  amount  of  high  Q  data.  The  low  amount  of  data  is  reflected  in  the 
high  Q  RSPs.  TM2  also  has  low  high  Q  RSP  values,  and  in  a  switch  from  the  previous 
lower  Q  maneuvers,  TM2  now  displays  worse  metrics  than  TM1  and  TM4,  both  of 
which  have  better  high  Q  RSPs.  The  better  RSPs  would  suggest  a  link  between  the 
high  Q  metrics  and  performance  of  the  maneuver  at  high  Q. 

Figure  39  shows  the  model  results  for  COM2.  Note  that  TM3  is  omitted.  The 
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Figure  36.  Model  results  COMla  Cl 


Figure  37.  Model  results  COMlb  Cl 
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Figure  38.  Model  results  COMlc  Cl 

omission  is  due  to  the  highly  inaccurate  curve  that  TM3  produced  which  distracted 
from  the  rest  of  the  results.  Also  consider  that  curve  TM4  is  essentially  showing  the 
fit  of  the  model  to  the  input  data,  rather  than  prediction  data.  TM1  and  TM2  are 
both  accurate  throughout  the  maneuver;  however,  TM2  comes  through  as  being  more 
accurate  in  the  metrics.  This  is  surprising  as  the  only  RSP  in  which  TM2  is  better 
than  TM1  is  percent  of  points  on  the  boundary.  The  success  of  TM2  would  suggest 
that  for  at  least  Cl,  hitting  the  boundaries  is  very  important. 

For  the  positive  Q  exponential  maneuver,  COM3,  shown  in  Figure  40,  all  of  the 
TMs  provided  good  results  through  the  course  of  the  maneuver.  However,  this  is  the 
first  instance  where  there  is  a  display  of  error  at  the  edge  of  the  regressor  space.  This 
error  by  TM3  occurs  because  that  portion  of  the  regressor  space  tested  by  COM3 
is  outside  the  bounds  of  TM3.  This  error  will  be  discussed  in  full  further  in  the 
document.  Fitting  with  the  pattern  from  the  previous  maneuvers,  TM2  returned  the 
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Figure  39.  Model  results  COM2  Cl 

best  values  while  TM1  returned  the  second  best  values  as  shown  in  Table  8.  These 
results  do  not  immediately  point  to  any  link  between  the  metrics  and  RSPs. 

The  results  for  COM4  shown  in  Figure  41  are  similar  to  those  from  the  previous 
maneuver.  The  models  all  line  up  well  in  the  linear  portion,  but  there  is  a  small  spread 
in  the  results  in  the  unstable  regime.  Since  this  maneuver  spends  a  longer  time  in 
the  unsteady  regime,  the  metrics  are  lower  overall  than  the  previous  maneuver.  From 
the  results  it  can  be  seen  that  TM3  is  close  but  not  accurate  in  the  high  alpha  low 
Q  portion  of  the  maneuver,  while  every  other  result  is  fairly  accurate.  TM3  also  had 
the  lowest  standard  deviation  value  for  the  low  Q  metric  cells.  This  could  indicate 
a  focus  point  of  not  having  a  well  spread  data.  Looking  at  the  tabular  results,  TM4 
has  better  metrics  than  TM2  for  the  most  accurate  maneuver. 

For  the  several  step  pitch  and  hold  maneuver,  COM5,  all  the  models  seem  to 
predict  the  zero  Q  values  to  a  relatively  high  degree  of  accuracy.  The  accuracy  of  the 
models  is  reflected  in  the  metrics;  however,  TM1  and  TM2  beat  the  other  maneuvers 
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Figure  40.  Model  results  COM3  Cl 


Figure  41.  Model  results  COM4  Cl 
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for  the  best  values  and  TM2  returns  better  values  than  TM1.  It  is  worth  noting  that 
TM3  loses  some  value  at  the  end  of  the  maneuver,  where  AoA  is  20  degrees,  likely 
due  the  fact  that  maximum  AoA  for  that  training  maneuver  is  around  18  degrees. 
Despite  the  large  difference  in  zero  Q  RSP  values,  there  does  not  seem  to  be  a  link 
to  the  metrics. 


Figure  42.  Model  results  COM5  Cl 


4.3.2  Coefficient  of  Drag 

The  results  for  Cd  can  be  seen  in  Table  9.  On  the  whole,  the  results  for  Cd  were 
less  accurate  than  the  results  seen  for  Cl-  There  is  also  a  larger  distinction  between 
the  different  training  maneuvers  than  was  seen  for  the  previous  coefficient. 

Looking  at  the  results  displayed  on  Figure  43,  it  is  easily  seen  that  TM1  and  TM4 
best  represent  the  data.  It  would  seem  that  now  that  the  data  is  not  as  well  organized 
as  lift  was,  the  ‘bad’  models  are  now  having  a  hard  time  accurately  predicting  the 
output.  TM4  resulted  in  the  largest  R2  and  lowest  NRMSD  as  shown  in  Table  9. 


Table  9.  R 2  and  NRMSD  for  all  Cd  models  vs  a  moving  average 


COM 

TM1 

TM2 

TM3 

TM4 

R2 

NRMSD 

R2 

NRMSD 

R2 

NRMSD 

R2 

NRMSD 

COMla 

0.9869 

0.0349 

0.9152 

0.1133 

0.8671 

0.0993 

0.9905 

0.0267 

COMlb 

0.9884 

0.0343 

0.8975 

0.1413 

0.9239 

0.0950 

0.9861 

0.0351 

COMlc 

0.9742 

0.0417 

0.9533 

0.0633 

0.7549 

0.2264 

0.9750 

0.0398 

COM2 

0.9809 

0.0430 

0.9115 

0.1146 

0.6048 

0.5714 

- 

- 

COM3 

0.7732 

0.0787 

0.8569 

0.0866 

0.7705 

0.0973 

0.9629 

0.0324 

COM4 

0.9668 

0.0621 

0.8475 

0.1567 

0.5973 

0.8259 

0.9836 

0.0414 

COM5 

0.9850 

0.0410 

0.9103 

0.1313 

0.9179 

0.1056 

0.9825 

0.0429 

TM4  has  the  largest  values  in  RSP1  and  RSP2  which  could  suggest  these  factors  are 
important  for  Cd  models. 

At  the  low-Q  AoA  sweep,  it  was  expected  that  TM3  would  still  do  well,  but  that 
proved  to  not  be  the  case  as  seen  in  Figure  44.  Both  TM3  and  TM2  seem  to  have 
a  problem  predicting  the  linear  stage,  whereas  in  the  unsteady  stage  they  are  much 
closer.  The  inability  of  these  two  maneuvers  to  predict  the  two  different  stages  could 
be  due  to  lack  of  a  variety  in  the  training  maneuvers  themselves.  From  the  metrics 
in  Table  9,  TM1  returns  better  values  than  TM4  in  both  R 2  and  NRMSD.  TM1  has 
better  low  Q  RSP  values,  and  since  this  is  low  Q  maneuver,  it  suggests  that  the  low 
Q  values  have  an  impact  on  the  metrics. 

As  expected  for  the  high  Q  sweep,  TM3  does  not  produce  very  good  results  as 
seen  in  Figure  45.  The  inaccuracy  makes  the  model  pretty  much  useless  in  these 
conditions.  TM1  and  TM4  again  display  good  values  throughout  the  majority  of  the 
maneuver,  only  beginning  to  separate  from  the  comparison  coefficient  at  the  high- 
alpha  values.  This  could  be  due  to  there  being  no  input  data  at  the  max  alpha,  at 
a  positive  Q  value.  In  the  training  maneuvers  at  the  max  alpha  the  maneuvers  are 
generally  at  a  near-zero  or  negative  Q  to  resulting  in  a  decrease  in  AoA.  TM4  returns 
the  best  values  for  this  comparison  maneuver.  TM4  also  has  the  best  high  Q  RSPs. 
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Figure  43.  Model  results  COMla  Cd 


Figure  44.  Model  results  COMlb  Cd 
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Figure  45.  Model  results  COMlc  Cd 

For  the  COM2  maneuver  shown  in  Figure  46,  TM3  again  was  so  inaccurate  that 
the  results  are  not  shown.  These  inaccuracies  likely  due  to  the  high  Q  nature  of 
this  maneuver  that  TM3  simply  does  not  have  the  data  for.  The  inaccuracy  of  TM3 
should  help  to  solidify  the  importance  of  covering  the  boundaries  of  the  regressor 
space.  TM4  also  is  not  taken  into  account  for  this  maneuver,  although  it  is  displayed 
on  the  figure.  From  Figure  46  and  Table  9  it  is  shown  that  TM1  is  the  most  accurate. 
TM1  does  have  better  high  Q  values  than  TM2  and  TM3,  but  it  is  hard  to  suggest 
that  these  factors  alone  are  why  it  does  better  for  this  maneuver. 

In  the  COM3  maneuver  shown  in  Figure  47,  there  is  a  large  spread  in  the  results. 
TM2  and  TM3  have  problems  in  the  low  alpha  and  low  Q  sections  of  the  maneuver. 
TM2  continues  to  display  inaccuracy  at  the  high  alpha  and  high  Q  portion  of  the  ma¬ 
neuver.  TM1  and  TM4  display  good  results  throughout  the  course  of  the  maneuver. 
By  the  tabular  results,  TM1  has  a  smaller  NRMSD  but  TM4  has  a  larger  R 2  value. 
The  results  are  very  close  in  magnitude.  TM1  and  TM4  have  better  low  Q  and  high 
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Figure  46.  Model  results  COM2  Cd 


Q  RSP  values  than  TM2  and  TM3,  but  TM3  still  manages  to  do  well  for  a  bit  of  the 
maneuver  at  high  Q.  There  is  no  clear  link  between  the  RSPs  and  metrics  for  this 
maneuver. 

For  the  COM4  maneuver  shown  in  Figure  48,  at  the  high  alpha,  low  Q  section 
of  the  maneuver  both  TM2  and  TM3  were  inaccurate.  TM3  retains  some  precision 
during  the  ramp  up  Q;  however,  it  again  becomes  inaccurate  at  the  higher  Q  portion 
of  the  maneuver.  TM1  and  TM4  both  retain  a  good  degree  of  accuracy  throughout 
the  course  of  the  maneuver.  The  only  lack  of  precision  being  the  very  minor  changes 
in  drag  at  low  angles  of  attack.  According  to  Table  9,  TM4  has  the  best  values  for 
R 2  and  NRMSD.  As  with  the  previous  maneuver,  it  is  hard  to  base  conclusions  on 
the  RSPs  from  this  result. 

For  the  last  comparison  maneuver,  which  is  the  step  maneuver  COM5,  the  results 
are  similar  to  the  previous  maneuvers.  TM2  and  TM3  have  trouble  at  the  lower 
alpha  values.  TM2  continues  to  stay  inaccurate  while  TM3  is  relatively  accurate  at 
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Figure  47.  Model  results  COM3  Cd 


Figure  48.  Model  results  COM4  Cd 
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15  degrees  AoA.  TM1  and  TM4  show  relative  accuracy  throughout  the  duration  of 
the  maneuver.  Table  9  shows  that  TM1  was  the  most  accurate  maneuver.  TM1  does 


have  the  best  zero  Q  RSP  results  for  percent  and  per  cell  but  TM4  has  a  better 
standard  deviation  number. 


Figure  49.  Model  results  COM5  Cn 


4.3.3  Moment  Coefficient 

The  results  for  the  moment  coefficient  models  can  be  seen  in  Table  10.  The 
prediction  of  results  for  the  moment  coefficient  fluctuated  for  the  models.  Overall, 
the  values  for  R 2  and  NRMSD  are  lower  than  the  values  seen  for  Cjj  and  Cl-  TM1 
and  TM4  have  the  best  RSP1  and  RSP 2  values. 

For  the  first  alpha  sweep,  COMla,  TM3  was  wildly  off  course  for  most  of  the 
maneuver,  whereas  TM1,  TM2  and  TM4  were  relatively  accurate.  TM1  is  a  bit 
inaccurate  at  the  higher  AoAs.  Per  the  metrics,  TM1  has  the  largest  R 2  while  TM4 
has  the  lowest  NRMSD. 
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Table  10.  R2  and  NRMSD  for  all  Cm  models  vs  a  moving  average 


COM 

TM1 

TM2 

TM3 

TM4 

R2 

NRMSD 

R2 

NRMSD 

R2 

NRMSD 

R2 

NRMSD 

COMla 

0.9487 

0.0556 

0.9096 

0.0624 

0.4612 

0.4436 

0.9434 

0.0504 

COMlb 

0.9366 

0.0871 

0.8642 

0.1073 

0.6078 

0.5391 

0.9182 

0.0902 

COMlc 

0.7634 

0.1351 

0.9106 

0.0854 

0.4263 

1.4006 

0.8825 

0.0879 

COM2 

0.9138 

0.0762 

0.8862 

0.0716 

0.5020 

0.2786 

- 

- 

COM3 

0.4243 

0.2198 

0.9252 

0.0589 

0.4668 

2.1220 

0.8438 

0.0641 

COM4 

0.8286 

0.1194 

0.8679 

0.0880 

0.5548 

1.0631 

0.9050 

0.0777 

COM5 

0.9118 

0.1015 

0.8960 

0.0954 

0.6196 

0.5177 

0.8851 

0.1108 

Figure  50.  Model  results  COMla  Cm 
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For  the  next  alpha  sweep,  COMlb,  TM3  is  again  not  very  accurate.  TM1  and 
TM4  show  the  same  trend  as  seen  in  the  previous  maneuver,  and  TM2  also  seems  to 
follow  the  same  line  as  before.  According  to  Table  10,  for  this  maneuver  TM1  has 
the  best  metrics.  As  seen  in  the  previous  coefficient,  TM1  has  the  best  metrics  and 
also  has  the  better  low  Q  RSPs.  Also  note  the  sharp  decrease  in  moment  coefficient. 
This  would  appear  to  correspond  to  separation.  While  there  is  a  discussion  to  be  had 
on  the  ability  of  CFD  to  accurately  predict  separation,  the  MVP  models  are  not  able 
to  capture  the  phenomena  that  is  predicted  by  the  CFD. 


Figure  51.  Model  results  COMlb  Cm 

For  the  final  alpha  sweep,  COMlc,  it  is  not  surprising  that  TM3  is  even  more 
inaccurate,  as  it  was  expected  to  be  for  the  high  Q  maneuvers.  Because  of  TM3s  in¬ 
accuracy  it  is  removed  from  Figure  52  as  it  was  distracting  from  the  other  maneuvers. 
As  it  can  be  seen,  both  TM4  and  TM1  have  a  sharp  rise  at  the  high  AoA  values.  Per 
Table  10,  TM2  has  the  best  output  metrics.  Much  of  the  inaccuracy  for  TM1  and 
TM4  seems  to  happen  in  the  last  0.1  seconds  of  the  maneuver;  without  this  section 
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TM1  and  TM4  would  likely  result  in  a  better  prediction.  Between  TM1  and  TM4, 
TM4  presents  the  better  metrics  and  has  better  RSP  values  at  high  Q. 


Figure  52.  Model  results  COMlc  Cm 

For  COM2,  TM4  should  not  be  considered  for  comparison  although  it  is  displayed 
in  Figure  53.  It  can  be  seen  that  TM1  and  TM2  both  overestimate  some  peaks  in  the 
unsteady  regime  while  overestimating  oscillations  in  the  linear  regime.  According  to 
Table  10,  TM1  is  the  most  accurate  maneuver.  TM1  does  have  better  high  Q  values 
than  TM2  and  TM3,  but  it  is  hard  to  suggest  that  these  factors  alone  are  why  it  does 
better  for  this  maneuver. 

For  COM3,  the  exponential  maneuver  shown  in  Figure  54,  TM3  is  not  shown 
due  to  gross  inaccuracy.  TM1  and  TM4  seem  to  capture  the  linear  regime  better; 
however,  TM2  seems  to  capture  the  non-linear  regime  better.  In  the  metrics,  TM4 
is  the  most  accurate  maneuver.  There  is  no  obvious  link  between  RSPs  and  output 
metrics  from  this  maneuver. 

For  the  negative  exponential  maneuver,  COM4,  shown  in  Figure  55,  TM3  is  accu- 
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Figure  53.  Model  results  COM2  Cm 


Figure  54.  Model  results  COM3  Cm 
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rate  in  the  low  Q  and  high  alpha  region.  However,  once  the  Q  value  increases,  TM3 
is  inaccurate.  TM1  and  TM4  are  again  closely  linked,  while  TM2  is  on  track  but 
seems  to  follow  the  trend  of  the  data  less  closely.  In  Table  10,  TM4  has  the  largest 
R2  while  NRMSD  has  the  lowest  NRMSD.  Once  again,  this  maneuver  provides  no 
obvious  insight  into  the  link  between  RSPs  and  the  metrics. 


Figure  55.  Model  results  COM4  Cm 

The  pitch  and  hold  maneuver  shown  in  Figure  56  does  not  display  TM3.  The 
results  of  TM3  were  so  inaccurate  that  it  was  difficult  to  distinguish  between  the 
other  results.  With  TM3  removed,  it  can  be  seen  that  TM4  and  TM1  do  well  at 
the  low  angles  of  attack,  have  trouble  around  10  degrees,  but  then  are  strangely  are 
better  at  15  and  20  degrees.  TM2  is  inaccurate  for  the  low  AoAs  but  does  well  at 
larger  angles.  Referencing  Table  10,  TM1  has  the  best  metrics  for  this  maneuver,  and 
as  before,  has  the  best  zero  Q  RSP  values. 
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Figure  56.  Model  results  COM5  Cm 


4.4  Statistical  Analysis  of  Initially  Proposed  RSPs 

The  results  of  this  section  and  the  results  of  Section  4.5.1  are  colored  using  the 
same  scale;  however,  each  coefficient  is  colored  using  a  different  scale.  The  differ¬ 
ent  coloring  is  due  to  the  large  discrepancy  between  the  correlation  values  for  each 
coefficient.  The  discrepancy  is  likely  due  to  the  scale  of  differences  between  R 2  and 
NRMSD  values  when  examining  coefficients.  For  Cl,  most  of  the  models  returned  a 
R 2  value  around  0.99,  while  for  Cm  a  high  R2  value  was  approximately  0.90,  while  a 
low  value  was  approximately  0.5.  The  correlation  values  are  then  smaller  for  Cl  as 
opposed  to  Cm  as  the  change  in  each  RSP  reflects  only  a  minor  at  best  difference  in 
R2  and  NRMSD  values.  Table  11  shows  the  coloring  for  the  correlation  values  based 
on  coefficient. 

Table  12  shows  the  colored  correlation  value  chart.  The  most  obvious  conclusion  is 
that  the  RSPs  are  not  consistently  correlated.  Despite  the  differing  ranges  to  provide 
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Table  11.  Correlation  value  coloring 


Color 

Correlation  Value  (Absolute  Value) 

cL 

Cd 

CM 

red 

x  <  0.5 

x  <  0.4 

x  <  0.6 

yellow 

0.5  <  x  <  0.6 

0.4  <  x  <  0.55 

0.6  <  x  <  0.8 

green 

x  >  0.6 

x  >  0.55 

x  >  0.8 

adequate  contrast,  the  correlation  values  are  not  consistent  across  coefficients.  That 
is,  a  good  RSP  for  one  coefficient  does  not  necessarily  lead  to  being  a  good  RSP  for 
another  coefficient.  Another  condition  of  note  is  that  the  standard  deviation  values 
seem  to  have  one  of  the  largest  impacts  across  most  the  coefficients.  Looking  at  RSP3 
and  RSP4,  it  appears  that  Q  is  a  larger  factor  than  AoA.  Therefore,  for  the  regressor 
space  subsets  (low  Q  and  high  Q)  it  may  be  beneficial  to  split  the  standard  deviation 
calculation  into  AoA  and  Q  standard  deviation  as  is  done  with  the  whole  regressor 
space  in  previous  RSPs. 


Table  12.  Initial  RSP  correlations  for  all  comparison  maneuvers 


RSP 

Description 

R2 

Cl 

NRMSD 

CD 

R2  NRMSD 

CM 

R2  NRMSD 

RSPi 

%  Whole 

0.3872 

-0.3738 

0.6027 

-0.4702 

0.4762 

-0.4220 

rsp2 

%  Boundary 

0.5279 

-0.5248 

0.5370 

-0.4268 

0.7129 

-0.5682 

rsp3 

a  std 

0.2821 

-0.2630 

0.4636 

-0.3824 

0.2873 

-0.3037 

rsp4 

Q  std 

0.4913 

-0.4694 

0.5860 

-0.5006 

0.5471 

-0.5230 

rsp5 

0Q  % 

0.3739 

-0.3562 

0.5742 

-0.4619 

0.4271 

-0.4045 

rsp6 

0Q  per  cell 

-0.3010 

0.3044 

-0.0785 

0.0906 

-0.4105 

0.3144 

rsp7 

0Q  std 

0.4713 

-0.4625 

0.6321 

-0.4902 

0.6168 

-0.5123 

rsp8 

LQ  % 

0.3135 

-0.2957 

0.5190 

-0.4180 

0.3435 

-0.3396 

rsp9 

LQ  per  cell 

-0.4685 

0.4655 

-0.2803 

0.2598 

-0.6040 

0.4934 

RSPi0 

LQ  std 

0.4519 

-0.4386 

0.6445 

-0.5063 

0.5641 

-0.4906 

RSPU 

HQ  % 

0.3966 

-0.3852 

0.6065 

-0.4685 

0.5025 

-0.4330 

RSPu 

HQ  per  cell 

0.2440 

-0.2334 

0.4931 

-0.3707 

0.3020 

-0.2710 

RSP13 

HQ  std 

0.5960 

-0.5828 

0.7018 

-0.5661 

0.7524 

-0.6414 

The  most  surprising  fact  is  that  none  of  the  values  for  Cl  or  Cm  are  very  strong. 
While  there  are  some  values  that  are  larger  than  others,  when  the  scale  is  held  con- 
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sistent  between  the  two  sections  (Section  4.4  and  4.5.1),  there  are  no  green  cells  for 
those  two  coefficients.  Cl  is  easy  to  predict,  causing  the  changes  in  RSPs  for  the 
different  training  maneuvers  to  have  little  effect  on  metrics.  Cm  is  difficult  for  all  the 
models  to  accurately  predict.  This  leads  to  the  impression  for  reformulation  of  the 
RSPs. 

The  first  reformulation  deals  with  how  these  standard  deviations  are  calculated. 
For  the  initial  RSPs,  the  distribution  is  filled  with  a  count  of  which  cells  are  hit.  A 
different  approach  is  to  use  the  total  number  of  points  in  a  cell,  rather  than  binary  1 
or  0  method. 

Since  measuring  the  shape  of  the  distribution  seems  to  be  important,  additional 
means  to  measure  the  distribution  might  provide  more  favorable  results. 

4.5  Proposed  RSPs 

After  the  analysis  conducted  in  Section  4.4,  new  RSPs  were  created  and  existing 
RSPs  were  modified.  There  were  also  several  changes  to  the  implementation  of  the 
programming,  which  in  some  cases  caused  a  small  difference  in  RSP  values  for  RSPs 
that  were  otherwise  not  changed. 

RSP\  is  unchanged  from  Equation  (3.4.1)  and  RSP-2  is  unchanged  from  Equation 
(3.4.2).  RSP3  no  longer  uses  the  maximum  value  as  a  standardization  and  takes  into 
account  how  many  times  a  value  is  hit.  This  is  in  order  to  conform  more  with  a 
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standard  statistical  cumulative  distribution  function.  The  revised  RSPs  is  shown  in 
Equation  (4.5.1). 


^  #  of  hits  at  first  AoA  ^ 

RSPs  =  std 

ff  of  hits  at  second  AoA 

4-  #  of  AoA  columns 

y  #  of  hits  at  last  AoA  y 

RSP4  makes  the  same  changes  in  regards  to  Q  as  RSPs  did  in  regards  to  AoA. 
The  revised  RSP4  is  given  in  Equation  (4.5.2). 


RSPa  =  std 


^  #  of  hits  at  first  Q  ^ 
ff  of  hits  at  second  Q 


[  \  #  °f  hits  at  last  Q  J 


4-  ft-  of  Q  rows 


(4.5.2) 


RSPs  is  unchanged  from  how  it  is  presented  in  Equation  (3.4.5)  and  RSPq  is 
unchanged  from  how  it  is  presented  in  Equation  (3.4.6). 

RSP7,  as  with  the  previous  standard  deviation  RSPs,  divides  by  the  number  of 
AoA  columns  as  shown  in  Equation  (4.5.3). 


RSP7  =  std 


cell2,  cells,  ...last  cell 
AoA  columns 


@0  Q 


(4.5.3) 


RSPs  is  unchanged  from  how  it  is  presented  in  Equation  (3.4.8),  and  RSPg  is 
unchanged  from  its  presentation  in  Equation  (3.4.9).  However  i?S'P10,  is  where  the 
RSPs  begin  to  lose  numbering  from  the  previous  iteration  of  RSPs.  RSP\o  is  now 
broken  up  into  two  RSPs,  which  collapse  the  low  Q  data  into  AoA  and  Q  cumulative 
distribution  functions  much  like  done  for  the  new  RSPs  and  RSP4  above.  The  new 
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RSP10  is  given  by  Equation  (4.5.4),  and  the  new  RSPn  is  given  by  Equation  (4.5.5). 


RSPW  =  std 


^  ff  of  hits  at  first  a  ^ 

ff  of  hits  at  second  a 

4-  jf  of  a  columns 

y  jf  of  hits  at  last  a  j 

(4.5.4) 


-I  @  low  ( not  0)  Q 


RSPn  =  std 


(  ff  of  hits  at  first  Q  ^ 

ff  of  hits  at  second  Q 

4-  //  of  Q  rows 

y  ff  of  hits  at  last  Q  j 

(4.5.5) 


-I  @  low  ( not  0)  Q 


The  high  Q  RSPs  are  set  up  the  same  as  the  new  low  Q  RSPs;  however,  the 
numbering  from  the  previous  iteration  of  RSPs  no  longer  lines  up.  The  new  RSPn  is 
equivalent  to  Equation  (3.4.11)  and  the  new  RSP\3  is  equivalent  to  Equation  (3.4.12). 
As  with  low  Q,  the  standard  deviation  is  now  split  into  two  different  standard  devi¬ 
ations,  given  in  Equations  (4.5.6)  and  (4.5.7). 


RSPia  =  std 


(  #  of  hits  at  first  a  ^ 

A  of  hits  at  second  a 

~  ff  of  a  columns 

y  ff  of  hits  at  last  a  j 

(4.5.6) 


J  @  high  Q 


RSPn  =  std 


#  of  hits  at  first  Q  ^ 

#  of  hits  at  second  Q 

4-  ff  of  Q  rows 

y  ff  of  hits  at  last  Q  j 

(4.5.7) 


-1  @  high  Q 


All  further  RSPs  presented  are  completely  unrelated  to  any  previous  RSPs  pre- 
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sented.  RSPig  is  an  attempt  to  measure  how  closely  and  completely  the  training 
maneuver  approaches  the  boundary.  For  the  RSP  calculations,  points  that  are  out¬ 
side  of  bounds  (±  75  deg/s  Q  and  -5  to  20  deg  a)  are  not  considered.  Since,  the 
maneuver  may  go  outside  the  bounds,  and  RSP16  aims  to  capture  how  many  points 
are  “wasted”.  RSPig  is  shown  in  Equation  (4.5.8). 


RSPW  = 


#  of  data  points  used 
T otal  #  of  points 


(4.5.8) 


The  next  two  RSPs  aim  to  quantify  how  efficiently  a  maneuver  captures  the  data. 
This  maneuver  takes  into  account  how  much  of  the  regressor  space  was  covered  and 
the  length  of  the  maneuver.  RSPn  is  given  by  Equation  (4.5.9)  and  RSPig  is  given 
by  Equation  (4.5.10).  RSPig  is  an  attempt  to  normalize  the  values  depending  on  the 
time  step  used  during  simulation  of  the  TMs. 


RSPn  = 
RSPig  = 


RSPi 

Length  of  M aneuver 
RSPi  *  dt 

Length  of  Maneuver 


(4.5.9) 

(4.5.10) 


The  next  three  RSPs  are  similar  to  the  previous  two  RSPs,  but  these  three  look 
to  the  how  much  of  each  subset  is  covered  per  the  total  duration  (length)  of  the 
maneuver  as  shown  in  Equations  (4.5.11)  -  (4.5.13). 


RSPi9  = 
RSP2  0  = 
rsp21  = 


rsp5 

Length  of  M aneuver 

RSPg 

Length  of  Maneuver 

RSPn 

Length  of  M aneuver 


(4.5.11) 

(4.5.12) 

(4.5.13) 


The  remaining  RSPs  are  different  from  the  previous  RSPs  in  that  they  are  gener- 
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ally  considered  higher  order  descriptions  of  data.  Three  different  higher  order  statis¬ 
tics  will  be  looked  at;  kurtosis,  skewness  and  Pearson  product-moment  correlation 
coefficient. 

Kurtosis  is  a  measure  of  the  “peakness”  of  a  data  set.  A  normal  distribution 
returns  a  kurtosis  value  of  3.  Kurtosis  is  defined  in  Equation  (4.5.14)  and  the  equation 
used  to  calculate  kurtosis  from  a  data  set  is  shown  in  Equation  (4.5.15).  [22] 


E(x-ii)4 

<j4 

n 

\  'ZiXi-x)4 

1=1 


n  \  ^ 

^E(^-^)2) 

i= 1  J 


(4.5.14) 

(4.5.15) 


where  x  is  the  independent  variable,  x  is  the  sample  mean,  yU  is  the  population  mean, 
n  is  the  number  of  data  points,  and  a  is  the  standard  deviation. 

Skewness  is  a  measure  of  the  asymmetry  of  the  data  around  the  mean.  A  positive 
skewness  represents  data  that  is  spread  out  more  to  the  right  of  the  mean.  Skewness 
is  defined  by  Equation  (4.5.16)  and  is  calculated  in  Equation  (4.5.17). 


s  = 


s  = 


E(x  —  /i)3 

<j3 

n 

±  £(Xi-x)3 


3 


(4.5.16) 

(4.5.17) 


The  Pearson  product-moment  correlation  coefficient  (Pearson  r)  is  a  measure  of 
linearity  between  two  data  sets.  A  value  of  1  means  positively  correlated  and  a  value 


of  -1  means  negatively  correlated.  The  equation  to  calculate  the  Pearson  r  is  given 
in  Equation  (4.5.18).  [26] 


r  = 


n 


E((®» 


x){yi  —  y)) 


n 


x)2J2(yi-y)2 

i=l 


(4.5.18) 


Using  the  above  statistics,  the  data  was  analyzed  in  the  various  subsets  as  was  done 
previously.  The  data  is  organized  into  two  different  distribution  curves,  one  for  AoA 
and  one  for  Q.  Those  distributions  are  also  partitioned  into  four  different  intervals: 
whole,  0  Q,  low  (non-0)  Q  and  high  Q.  RSP22  and  RSP2 3  look  at  the  kurtosis  of  the 
entire  regressor  space  in  Equation  (4.5.19)  and  Equation  (4.5.20). 


RSP22  =  Kurtosis 


^  #  of  hits  at  first  AoA  ^ 
ff  of  hits  at  second  AoA 


RSP2  3  =  Kurtosis 


y  ff  of  hits  at  last  AoA  J 

^  ff  of  hits  at  first  Q  ^ 
ff  of  hits  at  second  Q 


y  ff  of  hits  at  last  Q  J 


(4.5.19) 


(4.5.20) 


The  next  two  RSPs  look  at  the  skewness  of  the  entire  regressor  space  and  are 
shown  in  Equation  (4.5.21)  and  Equation  (4.5.22). 


RSP24  =  Skewness 


^  #  of  hits  at  first  AoA  ^ 
ff  of  hits  at  second  AoA 


y  A  of  Kts  at  last  AoA  J 


(4.5.21) 
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RSPoz  =  Skewness 


ff  of  hits  at  first  Q 
ff  of  hits  at  second  Q 


(4.5.22) 


ff  of  hits  at  last  Q 


The  next  RSP  is  the  Pearson  r  between  AoA  and  Q  for  the  maneuver  and  is  shown 


by  Equation  (4.5.23). 


RSP26  =  r  = 


^((AoA-AoA^Q.-Q)) 

2=1 

l~n  ^  n 

lUAoA-AoA^UQi-Q)2 

2=1  2=1 


(4.5.23) 


The  next  set  of  RSPs  looks  at  the  kurtosis  and  skewness  at  0  Q.  The  kurtosis  RSP 
is  given  by  Equation  (4.5.24)  and  the  skewness  RSP  is  given  by  Equation  (4.5.25). 


RSP?  7  =  Kurtosis 


ff  of  hits  at  first  AoA 
#  of  hits  at  second  AoA 


(4.5.24) 


ft  of  hits  at  last  AoA 


RSP?r  =  Skewness 


ff  of  hits  at  first  AoA 
ff  of  hits  at  second  AoA 


(4.5.25) 


ff  of  hits  at  last  AoA 


The  next  set  of  equations  looks  at  the  kurtosis  and  skewness  for  Q  and  AoA  at 
low  but  not  0  Q.  Kurtosis  for  AoA  is  given  in  Equation  (4.5.26),  and  Q  is  given  in 


Equation  (4.5.27).  Skewness  for  AoA  is  given  in  Equation  (4.5.28),  and  Q  is  given  in 
Equation  (4.5.29). 


RSPo  q  =  Kurtosis 


RSP30  =  Kurtosis 


(  ft  of  hits  at  first  AoA  ^ 
ft  of  hits  at  second  AoA 


y  ft  of  hits  at  last  AoA  J 


^  ft  of  hits  at  first  Q  ^ 
ft  of  hits  at  second  Q 


@  low  ( not  0)  Q 


RSP31  =  Skewness 


y  ff  of  hits  at  last  Q  J 

^  #  of  hits  at  first  AoA 
ft  of  hits  at  second  AoA 


low  ( not  0)  Q 

\ 


RSP32  =  Skewness 


y  ft  of  hits  at  last  AoA  J 


(  ft  of  hits  at  first  Q  ^ 
ft  of  hits  at  second  Q 


@  low  ( not  0)  Q 


y  #  of  hits  at  last  Q  J 


@  low  ( not  0)  Q 


(4.5.26) 


(4.5.27) 


(4.5.28) 


(4.5.29) 


The  next  set  of  RSPs  measures  the  kurtosis  and  skewness  of  the  high  Q  data. 
Equation  (4.5.30)  measures  the  kurtosis  of  AoA  and  Equation  (4.5.31)  measures  the 


kurtosis  of  Q.  Equation  (4.5.32)  measures  the  skewness  of  AoA  and  Equation  (4.5.33) 
measures  the  skewness  of  Q. 


RSP33  =  Kurtosis 


RSP3  4  =  Kurtosis 


^  if  of  hits  at  first  AoA  ^ 
ff  of  hits  at  second  AoA 


y  if  of  hits  at  last  AoA  J 

^  ff  of  hits  at  first  Q  ^ 
if  of  hits  at  second  Q 


@  high  Q 


RSP35  =  Skewness 


y  if  of  hits  at  last  Q  J 

(  if  of  hits  at  first  AoA 
if  of  hits  at  second  AoA 


high  Q 

\ 


RSP3  K  =  Skewness 


y  if  of  hits  at  last  AoA  J 


(  if  of  hits  at  first  Q  ^ 
if  of  hits  at  second  Q 


)  high  Q 


y  #  of  hits  at  last  Q  J 


)  high  Q 


(4.5.30) 


(4.5.31) 


(4.5.32) 


(4.5.33) 


RSP36  concludes  the  RSPs  considered  for  this  research.  The  MATLAB®function 
created  for  the  calculation  of  these  RSPs  is  given  in  Appendix  A. 4  for  the  reader’s 
reference. 
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4.5.1  Proposed  RSP  Results 


Following  the  discussion  in  Section  4.5,  the  additional  RSPs  were  calculated  in 
two  groups.  First,  the  ‘first  order’  statistics  are  given  in  Table  13.  The  first  thing  to 
note  is  the  different  standard  deviation  values.  While  it  was  expected  for  the  values 
to  change  as  the  normalization  was  changed,  the  change  in  switching  to  total  hits 
instead  of  cells  has  changed  the  trend  of  the  data.  The  change  is  most  evident  in 
RSPs  3,  4  and  7.  In  Table  6,  TM1  has  the  largest  AoA  standard  deviation,  while  in 
Table  13  it  has  the  third  largest.  The  implications  of  these  different  results  may  be 
seen  later  on  with  the  correlation  results. 


Table  13.  First  set  of  revised  RSP  values 
RSP  Description  TM1  TM2  TM3  TM4 


RSP1  %  Whole 
RSP2  %  Boundary 


RSP3 

a  std 

RSP4 

Q  std 

RSP5 

OQ  % 

rsp6 

OQ  per  cell 

rsp7 

OQ  std 

RSPs 

LQ  % 

rsp9 

LQ  per  cell 

RSPW 

LQ  a  std 

RSPn 

LQ  Q  std 

rsp12 

HQ  % 

RSP13 

HQ  per  cell 

RSPU 

HQ  a  std 

RSP15 

HQ  Q  std 

rsp16 

Effectiveness 

rsp17 

Tirnel 

RSPis 

Time2 

RSP19 

Time  Eff  OQ 

RSP20 

Time  Eff  LQ 

RSP2 1 

Time  Eff  HQ 

3.6573 

0.8851 

3.4897 

11.556 

3.6625 

1.3754 

0.0861 

0.0098 

4.8000 

0.8000 

0.8720 

0.1280 

0.0178 

0.0057 

4.9560 

0.8120 

0.9473 

0.1311 

2.9953 

1.3198 

0.2047 

0.0024 

3.1866 

0.9115 

0.8153 

0.1874 

1.9957 

0.6853 

0.1285 

0.0142 

0.9964 

0.8082 

0.9143 

0.8851 

1.142e-05 

1.106e-05 

1.2000 

0.8000 

1.2390 

0.8120 

0.7966 

0.9115 

0.6880 

5.0024 

0.0000 

15.046 

4.3837 

4.9181 

0.2083 

0.0935 

0.8000 

4.4000 

0.9880 

0.7640 

0.0488 

0.0160 

1.3050 

4.2230 

1.3432 

0.8761 

4.4692 

4.4361 

1.1168 

0.0158 

0.4650 

5.2843 

0.3845 

1.1357 

1.1915 

2.6407 

0.1935 

0.1456 

1.0000 

0.9999 

0.2293 

1.0005 

2.866e-06 

1.251e-05 

0.2667 

0.8800 

0.4350 

0.8446 

0.1550 

1.0569 

Another  note  is  the  switch  from  one  standard  deviation  value  for  the  low  Q  and 
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high  Q  ranges  to  one  for  AoA  and  one  for  Q.  In  Table  13,  it  can  be  seen  that  there 
are  different  trends  in  AoA  and  Q. 

A  final  note  concerns  the  trends  seen  in  RSPs  19,  20  and  21,  these  RSPs  are  far 
different  than  the  percentage  RSPs  (RSPs  5,8  and  12)  that  form  the  basis.  The  result 
of  this  can  be  seen  in  the  magnitude  of  the  correlation  values. 

Table  14  shows  the  correlation  values  for  the  ‘first  order’  statistics.  The  positive 
outlook  is  now  we  can  see  green  start  to  cover  more  of  the  table.  Some  statistics  are 
starting  to  reveal  themselves  as  more  useful  than  others. 


Table  14.  First  set  of  revised  RSPs  -  correlation  values 


RSP 

Description 

R2 

cL 

NRMSD 

cD 

R2  NRMSD 

CM 

R2  NRMSD 

RSP1 

%  Whole 

0.3838 

-0.3648 

0.6001 

-0.4666 

0.4626 

-0.4166 

RSPo 

%  Boundary 

0.4889 

-0.4827 

0.4920 

-0.3862 

0.6644 

-0.5253 

RSP3 

a  std 

-0.2577 

0.2610 

-0.0015 

0.0502 

-0.3242 

0.2649 

RSP4 

Q  std 

-0.6034 

0.5936 

-0.4752 

0.4207 

-0.7563 

0.6386 

rsp5 

0Q  % 

0.3716 

-0.3503 

0.5721 

-0.4595 

0.4181 

-0.4009 

rsp6 

0Q  per  cell 

-0.3660 

0.3686 

-0.1534 

0.1538 

-0.4923 

0.3845 

rsp7 

0Q  std 

-0.6474 

0.6346 

-0.5691 

0.4901 

-0.8098 

0.6874 

RSP8 

LQ  % 

0.3117 

-0.2910 

0.5189 

-0.4166 

0.3370 

-0.3369 

rsp9 

LQ  per  cell 

-0.4728 

0.4707 

-0.2851 

0.2641 

-0.6141 

0.4985 

RSPo 

LQ  a  std 

-0.3760 

0.3740 

-0.1462 

0.1717 

-0.4591 

0.3912 

RSPU 

LQ  Q  std 

-0.6634 

0.6484 

-0.6405 

0.5384 

-0.8311 

0.7068 

rsp12 

HQ  % 

0.3930 

-0.3754 

0.6036 

-0.4648 

0.4872 

-0.4271 

rsp13 

HQ  per  cell 

0.2405 

-0.2241 

0.4903 

-0.3672 

0.2863 

-0.2653 

RSPU 

HQ  a  std 

0.2069 

-0.1911 

0.4622 

-0.3427 

0.2441 

-0.2297 

rsp15 

HQ  Q  std 

-0.4483 

0.4463 

-0.2409 

0.2357 

-0.5731 

0.4708 

RSPi6 

Effectiveness 

-0.2523 

0.2584 

-0.0067 

0.0397 

-0.3457 

0.2615 

rsp17 

Timel 

0.6652 

-0.6472 

0.6883 

-0.5745 

0.8218 

-0.7096 

RSPis 

Time2 

0.6652 

-0.6472 

0.6883 

-0.5745 

0.8218 

-0.7096 

rsp19 

Time  Eff  0Q 

-0.3520 

0.3310 

-0.5577 

0.4466 

-0.3932 

0.3802 

rsp20 

Time  Eff  LQ 

-0.6577 

0.6435 

-0.5909 

0.5094 

-0.8153 

0.6983 

rsp21 

Time  Eff  HQ 

-0.2936 

0.2731 

-0.5024 

0.4032 

-0.3128 

0.3175 

The  first  statistic  in  this  group  is  the  time  statistics.  For  this  research  a  constant 
time  step,  At  was  used,  hence,  RSP\7  and  RSPi$  have  the  same  correlation  values. 
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It  is  interesting  that  these  RSPs  would  be  more  significant  than  RSP\ ,  which  they 
are  based  on. 

The  zero  Q  standard  deviation  has  become  much  more  useful  due  to  the  restruc¬ 
turing  of  calculation.  The  thought  behind  the  zero  Q  data  was  that  it  might  be 
important  how  long  a  maneuver  spends  in  each  cell,  hence  the  creation  of  RSPq. 
However,  by  calculating  the  standard  deviation  based  on  the  amount  of  the  points  at 
that  AoA  as  opposed  to  by  a  cell  count,  that  ‘length  of  hover’  is  now  incorporated 
into  the  statistic  as  well  as  the  spread  across  the  AoA  regime. 

The  separation  of  AoA  and  Q  for  the  standard  deviation  calculations  and  the 
overall  switch  to  number  of  total  hits  versus  number  of  cell  hits,  seems  to  have  done 
well  for  the  low  Q  metrics.  With  the  revised  setup,  low-Q  standard  deviation  of  Q 
has  been  shown  to  be  a  significant  factor,  whereas  low  Q  AoA  standard  deviation  as 
been  shown  as  insignificant.  This  implies  Q  is  a  greater  consideration  than  AoA  for 
the  investigated  conditions 

While  the  next  system  seemed  to  work  well  for  the  low  Q  values,  the  high  Q  values 
do  not  seem  to  have  improved.  Neither  the  AoA  or  the  Q  standard  deviation  has  an 
improved  correlation  to  the  old  system  of  calculation. 

The  second  set  of  RSPs  is  shown  in  Table  15.  One  interesting  statistic  is  that 
RSP-26  is  nearly  zero  for  all  cases.  Another  interesting,  or  at  least  promising  trend, 
is  that  for  each  subsection  there  is  a  difference  in  the  values.  This  should  insure  that 
there  are  some  different  correlations. 

The  correlation  values  for  the  second  set  of  RSPs  is  shown  in  Table  16.  The  first 
thing  to  notice  is  the  chart  has  much  more  yellow  and  green  correlations  than  red. 
The  change  in  coloring  shows  the  revised  RSPs  have  been  determined  to  be  more 
useful  when  compared  to  the  previous  set  of  RSPs.  Kurtosis  for  AoA,  LQ  skewness 


93 


Table  15.  Second  set  of  final  RSP  values 


RSP 

Description 

TM1 

TM2 

TM3 

TM4 

rsp22 

Kurtosis  a 

8.8673 

45.752 

44.091 

12.483 

psp23 

Kurtosis  Q 

27.190 

7.1970 

108.31 

29.246 

rsp24 

Skewness  a 

2.0055 

5.4952 

4.9727 

2.6874 

RSP25 

Skewness  Q 

3.0550 

2.0665 

8.4767 

3.7456 

RSP26 

Pearson 

0.0000 

0.0000 

0.0000 

0.0458 

rsp27 

Kurtosis  OQ 

56.051 

123.01 

197.99 

34.044 

RSP2s 

Skewness  OQ 

6.6218 

11.046 

13.700 

5.5731 

rsp29 

Kurtosis  LQ  a 

9.0705 

70.075 

48.588 

17.429 

RSPso 

Kurtosis  LQ  Q 

39.463 

1.2336 

68.916 

3.4526 

RSP3 1 

Skewness  LQ  a 

2.3732 

7.8303 

5.6304 

3.7178 

RSP32 

Skewness  LQ  Q 

4.9504 

0.4401 

7.5377 

-0.2166 

RSP33 

Kurtosis  HQ  a 

2.4077 

2.8856 

1.6763 

4.1160 

RSp34 

Kurtosis  HQ  Q 

24.574 

5.7243 

120.328 

21.999 

RSP35 

Skewness  HQ  a 

0.3376 

0.4104 

-0.4953 

1.0433 

RSP36 

Skewness  HQ  Q 

2.9407 

1.7515 

8.0510 

3.1100 

of  AoA,  high  Q  kurtosis  of  AoA  and  high  Q  skewness  of  Q  are  all  determined  to  be 
significant  for  each  coefficient. 

The  task  now  remains  to  reduce  the  list  of  36  RSPs  to  an  easier  to  manage  list 
of  RSPs  from  which  to  motivate  a  maneuver.  From  Section  4.3  it  was  determined 
the  zero  Q,  low  Q  and  high  Q  results  trended  towards  more  accurate  maneuvers  in 
those  areas.  Thus,  it  tracks  that  at  least  one  RSP  to  be  maximized  for  the  ‘ideal’ 
maneuver.  RSPn  has  some  of  the  highest  correlation  values  across  the  board  and 
is  a  measure  of  how  fast  that  space  is  covered.  Therefore,  it  is  fitting  RSPn  is  the 
primary  RSP  to  use  in  designing  a  maneuver. 

The  next  step  is  to  determine  which  low  Q  RSP  to  use.  RSP34  and  RSPn  return 
very  high  correlation  values,  and  seem  to  be  strong  measures  of  maneuver  effective¬ 
ness.  However,  RSPn  has  slightly  higher  correlation  values  for  most  of  the  metrics, 
and  is  therefore  selected  to  go  forward. 

For  the  high  Q  RSP  selection,  the  choice  is  equally  as  tough.  RSP3 3  and  RSP36 
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Table  16.  Second  set  of  revised  RSPs  -  correlation  values 


RSP 

Description 

CL 

R 2  NRMSD 

cD 

R2  NRMSD 

CM 

R2  NRMSD 

rsp22 

Kurtosis  a 

-0.6548 

0.6401 

-0.5826 

0.5065 

-0.8051 

0.6946 

RSP23 

Kurtosis  Q 

0.2123 

-0.2055 

0.3994 

-0.2758 

0.3126 

-0.2371 

RSP2A 

Skewness  a 

-0.5723 

0.5513 

-0.7033 

0.5715 

-0.6936 

0.6140 

rsp25 

Skewness  Q 

-0.5180 

0.4968 

-0.6794 

0.5476 

-0.6214 

0.5569 

RSP26 

Pearson 

-0.1751 

0.1569 

-0.4088 

0.3189 

-0.1676 

0.1923 

rsp27 

Kurtosis  OQ 

-0.5696 

0.5617 

-0.5189 

0.4261 

-0.7496 

0.6081 

RSP2s 

Skewness  OQ 

-0.1207 

0.1032 

-0.3461 

0.2717 

-0.0918 

0.1334 

RSP-29 

Kurtosis  LQ  a 

-0.5328 

0.5260 

-0.5038 

0.4054 

-0.7123 

0.5704 

RSP30 

Kurtosis  LQ  Q 

0.4829 

-0.4724 

0.5739 

-0.4429 

0.6448 

-0.5215 

RSP31 

Skewness  LQ  a 

-0.6650 

0.6500 

-0.6170 

0.5272 

-0.8250 

0.7069 

RSP32 

Skewness  LQ  Q 

0.5852 

-0.5700 

0.6665 

-0.5326 

0.7494 

-0.6286 

RSP33 

Kurtosis  HQ  a 

-0.6585 

0.6441 

-0.5930 

0.5113 

-0.8155 

0.6992 

RSP3  4 

Kurtosis  HQ  Q 

0.5959 

-0.5743 

0.6383 

-0.5480 

0.6930 

-0.6335 

RSP35 

Skewness  HQ  a 

0.5375 

-0.5160 

0.5792 

-0.5045 

0.6072 

-0.5703 

RSP36 

Skewness  HQ  Q 

0.6499 

-0.6345 

0.6655 

-0.5491 

0.8202 

-0.6943 

correlated  to  the  metrics  better  than  any  of  the  other  high  Q  RSPs.  RSP^q  has  larger 
correlation  values  than  RSP33,  so  it  was  chosen  to  be  optimized. 

The  next  step  is  to  look  at  which  zero  Q  RSP  to  optimize.  RSP7  returned  the 
highest  correlation  values  by  a  good  margin,  so  it  was  selected  to  be  optimized. 

Finally,  it  must  be  decided  which  RSP  will  be  most  important.  There  are  a  lot 
of  parameters  to  vary  and  it  seems  unlikely  there  is  one  maneuver  that  will  have  the 
best  values  for  all  these  parameters.  RSPn  has  very  high  correlation  values  and  is  a 
measure  of  the  entire  maneuver  it  is  the  most  important.  However,  since  low  Q  values 
are  more  difficult  to  predict,  RSPn  is  selected  as  the  next  parameter.  High  Q  values 
are  a  very  large  subset  of  the  total  regressor  space  so  RSP36  is  next.  Finally,  since 
zero  Q  data  is  the  smallest  subset  of  the  regressor  space,  it  is  the  least  important. 
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4.6  RSP  Motivated  Training  Maneuver 


Following  the  discussion  in  Section  4.5.1,  RSPs  17,  11,  36,  and  7  were  decided  to 
be  used  to  motivate  a  training  maneuver.  The  training  maneuver  was  first  attempted 
as  an  expansion  of  the  “chirp”  maneuver.  However,  it  was  found  to  be  difficult  to 
improve  the  RSPs  using  a  “chip”  type  maneuver.  Therefore,  a  training  maneuver 
similar  to  TM4  was  adopted.  TM4  is  the  training  maneuver  adopted  from  COM2 
in  order  to  provide  another  set  of  RSPs  for  the  statistical  analysis.  The  simularity 
can  be  seen  in  Figure  57.  In  order  to  increase  RSPn,  the  length  of  the  maneuver 
was  shortened.  While  TM5  looks  similar  to  TM4,  is  it  quite  different.  The  difference 
can  be  seen  most  clearly  in  Figure  58.  TM4  utilizes  changes  in  both  amplitude  and 
frequency  for  both  sinusoids,  whereas  TM4  has  constant  amplitude  and  frequency 
for  both  sinusoid  functions.  Lastly  note  the  large  spread  of  data  across  the  regressor 
space  shown  in  Figure  59.  It  was  also  quickly  determined  that  hireling  a  maneuver 
which  has  the  best  value  for  each  of  the  selected  RSPs  is  extremely  difficult. 

The  selected  RSPs  are  calculated  for  each  training  maneuver  and  shown  in  Table 
17.  The  table  also  shows  whether  the  RSP  should  be  maximized  or  minimized.  From 
the  table,  it  can  be  seen  that  TM5  ranks  the  best  RSPn,  the  third  best  RSPn,  the 
second  best  RSP36  and  the  second  best  RSP7.  If  the  RSPs  are  truely  and  indication 
of  future  training  maneuver  performance,  these  RSP  values  would  point  to  TM5 
producing  accurate  results. 


Table  17.  RSP  motivated  TM  and  RSP  values 


RSP 

Description 

Max  or  Min 

TM1 

TM2 

TM3 

TM4 

TM5 

RSP17 

Time 

Max 

0.9143 

0.8851 

0.2293 

1.0005 

1.2660 

RSPn 

LQ  Q  std 

Min 

0.2047 

0.0024 

1.1168 

0.0158 

0.0234 

RSP35 

Skewness  HQ  Q 

Min 

2.9407 

1.7515 

8.0510 

3.1100 

2.3560 

rsp7 

OQ  std 

Min 

0.0178 

0.0057 

0.0488 

0.0160 

0.0116 
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Figure  57.  Training  maneuver  5:  AoA  and  Q  vs  time 


AoA  (deg) 

Figure  58.  Training  maneuver  5:  Q  vs  AoA 
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Pitch  Rate  (Q) 


Pitch  Rate  (Q)  _ioo  _io 


AoA 


Figure  59.  Training  maneuver  5:  discretized  regressor  space 


4.6.1  RSP  Motivated  Maneuver  Results 

In  an  effort  to  reduce  presentation  of  redundant  information,  the  metrics  for  TM5 
are  not  presented  next  to  the  metric  information  for  the  other  training  maneuvers. 
Instead  the  data  will  be  presented  for  all  comparison  maneuvers  and  all  coefficients 
for  just  TM5.  The  chart  is  color-coded  podium  style,  with  gold  symbolizing  the  ‘best’ 
comparison  metrics  when  compared  to  the  other  four  models,  silver  second  and  bronze 
third.  This  data  is  represented  in  Table  18. 

In  another  effort  to  reduce  the  amount  of  redundant  information  presented,  only 
charts  deemed  important  to  the  understanding  of  TM5  model  performance  will  be 
presented  here.  However,  the  charts  for  each  coefficient  and  comparison  maneuver 
with  TM5  plotted  are  included  in  Appendix  B. 

As  can  be  seen  from  Table  18,  TM5  does  very  well  for  some  maneuvers,  accurate 
for  some  others  and  then  comparably  for  the  rest.  There  are  no  points  at  which  TM5 
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Table  18.  TM5  model  metrics 


COM 

R2 

CL 

NRMSD 

R2 

cD 

NRMSD 

CM 

R2  NRMSD 

COMla 

0.9875 

0.0292 

0.9878 

0.0301 

0.9182 

0.0592 

COMlb 

0.9843 

0.0325 

0.9785 

0.0435 

0.8674 

0.1145 

COMlc 

0.9873 

0.0264 

0.9640 

0.0456 

0.8088 

0.1050 

COM2 

0.9862 

0.0279 

0.9869 

0.0337 

0.9236 

0.0778 

COM3 

0.9950 

0.0185 

0.8927 

0.0516 

0.5881 

0.1123 

COM4 

0.9824 

0.0280 

0.9865 

0.0370 

0.9200 

0.0658 

COM5 

0.9784 

0.0354 

0.9755 

0.0512 

0.8371 

0.1364 

is  a  terrible  maneuver.  Looking  at  the  cells  that  are  not  colored  at  all,  we  see  that 
the  R2  values  are  large  and  NRMSD  values  are  small,  indicating  that  all  models  were 
extremely  accurate  for  those  maneuvers,  and  the  difference  among  all  the  models  is 
slight.  Just  because  TM5  did  not  get  the  best  values  in  those  cases  does  not  mean  it 
is  an  unsuitable  TM. 

The  results  for  COMlb  for  Cl  are  shown  in  Figure  60.  While  TM5  did  not  place 
in  the  top  three  results,  it  produces  highly  accurate  results  throughout  the  range 
AoAs. 

Figure  61  shows  the  high  accuracy  of  the  TM5  model.  While  most  of  the  other 
models  had  trouble  resolving  the  low  AoA  high  Q  aspect  of  COM4,  TM5  was  able  to 
retain  its  fidelity.  This  ability  seems  to  suggest  the  favorable  value  of  RSP- 36  leads  to 
good  high  Q  results. 

The  one  area  of  concern,  which  is  a  consistent  area  of  concern  for  all  the  models,  is 
COM3  for  Cm-  TM5  returns  a  R2  value  of  0.5881.  The  model  is  shown  in  Figure  62. 
The  reason  for  the  low  R2  value  is  the  high  Q,  high  alpha  portion  of  COM3.  While  the 
model  has  high  Q  and  high  alpha  data  points,  it  does  not  have  a  good  representation 
of  both  of  them.  Thinking  back  to  the  regressor  space  figures  in  Section  4.6,  there  is 
very  little,  or  no  data,  in  the  corners  of  the  regressor  space.  When  the  model  is  then 
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Figure  60.  Model  results  COMlb  Cl  with  TM5 


Figure  61.  Model  results  COM4  Cd  with  TM5 
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mapped  to  the  resulting  data  of  the  training  maneuver,  this  subset  of  data  is  outside 
the  bounds  of  the  training  maneuver,  and  the  polynomial  functions  fall  off. 


Figure  62.  Model  results  COM3  Cm  with  TM5 

At  this  point,  it  has  been  shown  using  statistical  measures  to  gain  greater  insight 
into  actions  of  a  training  maneuver,  the  maneuver  can  be  better  tailored  to  a  partic¬ 
ular  need.  However,  with  there  being  some  time  left  an  additional  maneuver  may  be 
able  to  solidify  or  expand  the  knowledge  of  these  RSPs  and  the  output  metrics.  From 
the  experience  gained  over  the  course  of  this  research,  the  author  has  decided  a  series 
of  actions  that  can  help  with  the  system  identification  process.  First,  the  results  from 
TM5  should  be  appended  to  the  statistical  analysis  to  better  determine  which  RSPs 
are  important.  Applying  the  results  of  TM5  to  the  data  set  and  recalculating  the 
correlation  values  gives  the  chart  shown  in  Table  19. 

Comparing  the  results,  most  of  the  important  RSPs  from  before  are  decidedly 
still  useful  to  predicting  TM  accuracy.  However,  the  values  of  the  correlations  have 
changed.  From  the  results,  RSPn  now  has  stronger  correlation  values  than  RSPn , 
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Table  19.  Correlation  values  with  TM5  included 


RSP 

Description 

cL 

R2  NRMSD 

CD 

R2  NRMSD 

CM 

R2  NRMSD 

RSP] 

%  Whole 

0.4002 

-0.3740 

0.6153 

-0.4891 

0.4223 

-0.4358 

rsp2 

%  Boundary 

0.4520 

-0.4498 

0.4240 

-0.3421 

0.6487 

-0.4841 

RSP3 

a  std 

-0.2534 

0.2572 

-0.0006 

0.0483 

-0.3217 

0.2619 

rsp4 

Q  std 

-0.6166 

0.6043 

-0.4994 

0.4466 

-0.7492 

0.6547 

rsp5 

0Q  % 

0.3917 

-0.3637 

0.5955 

-0.4846 

0.3756 

-0.4247 

rsp6 

0Q  per  cell 

-0.3532 

0.3574 

-0.1347 

0.1412 

-0.4980 

0.3716 

rsp7 

0Q  std 

-0.6595 

0.6431 

-0.5970 

0.5199 

-0.7927 

0.7025 

rsp8 

LQ  % 

0.3418 

-0.3139 

0.5499 

-0.4464 

0.3067 

-0.3715 

rsp9 

LQ  per  cell 

-0.4775 

0.4748 

-0.2917 

0.2738 

-0.6228 

0.5049 

RSP\o 

LQ  a  std 

-0.3948 

0.3901 

-0.1860 

0.2035 

-0.4614 

0.4141 

RSPn 

LQ  Q  std 

-0.6741 

0.6553 

-0.6636 

0.5660 

-0.8099 

0.7199 

RSP12 

HQ  % 

0.4113 

-0.3863 

0.6232 

-0.4919 

0.4539 

-0.4480 

RSPi3 

HQ  per  cell 

0.2858 

-0.2634 

0.5378 

-0.4154 

0.2952 

-0.3141 

RSPU 

HQ  a  std 

0.2627 

-0.2394 

0.5174 

-0.3976 

0.2603 

-0.2901 

RSP15 

HQ  Q  std 

-0.4595 

0.4561 

-0.2621 

0.2556 

-0.5813 

0.4849 

rsp16 

Effectiveness 

-0.2433 

0.2504 

0.0018 

0.0323 

-0.3562 

0.2532 

RSP17 

Timel 

0.6244 

-0.5995 

0.6867 

-0.5767 

0.7121 

-0.6700 

RSPis 

Time2 

0.6244 

-0.5995 

0.6867 

-0.5767 

0.7121 

-0.6700 

rsp19 

Time  Eff  0Q 

0.5737 

-0.5459 

0.6512 

-0.5586 

0.5956 

-0.6140 

RSP20 

Time  Eff  LQ 

0.5070 

-0.4789 

0.5922 

-0.5099 

0.4987 

-0.5430 

RSP2i 

Time  Eff  HQ 

0.6284 

-0.6060 

0.6793 

-0.5662 

0.7439 

-0.6746 

RSP-22 

Kurtosis  a 

-0.3828 

0.3553 

-0.5933 

0.4824 

-0.3648 

0.4149 

RSP23 

Kurtosis  Q 

-0.6688 

0.6508 

-0.6183 

0.5389 

-0.7915 

0.7124 

RSP24 

Skewness  a 

-0.3270 

0.2992 

-0.5355 

0.4343 

-0.2869 

0.3558 

RSP25 

Skewness  Q 

-0.6669 

0.6491 

-0.6062 

0.5332 

-0.7818 

0.7094 

rsp26 

Spearman 

0.2083 

-0.2021 

0.3790 

-0.2683 

0.3263 

-0.2311 

rsp27 

Kurtosis  0Q 

-0.5685 

0.5416 

-0.7127 

0.5876 

-0.6289 

0.6121 

rsp28 

Skewness  0Q 

-0.5106 

0.4830 

-0.6817 

0.5573 

-0.5488 

0.5513 

rsp29 

Kurtosis  LQ  a 

-0.2312 

0.2064 

-0.4709 

0.3753 

-0.1734 

0.2534 

RSP30 

Kurtosis  LQ  Q 

-0.5811 

0.5671 

-0.5622 

0.4675 

-0.7416 

0.6232 

RSP3 1 

Skewness  LQ  a 

-0.1844 

0.1601 

-0.4187 

0.3350 

-0.1050 

0.2028 

RSP32 

Skewness  LQ  Q 

-0.5473 

0.5340 

-0.5497 

0.4494 

-0.7111 

0.5884 

RSP33 

Kurtosis  HQ  a 

0.4410 

-0.4349 

0.4935 

-0.3908 

0.6148 

-0.4741 

rsp34 

Kurtosis  HQ  Q 

-0.6757 

0.6570 

-0.6419 

0.5554 

-0.8005 

0.7204 

RSP35 

Skewness  HQ  a 

0.5812 

-0.5668 

0.6424 

-0.5262 

0.7374 

-0.6233 

RSP36 

Skewness  HQ  Q 

-0.6699 

0.6519 

-0.6193 

0.5400 

-0.7923 

0.7135 
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and  shifting  from  previous  results,  RSPU  has  a  greater  correlation  than  RSP36. 
Therefore,  the  decided  new  order  of  RSP  importance  is  RSPn ,  RSPn,  RSP34  and 
RSP7. 

Second,  the  phenomena  discussed  above  were  observed  for  many  models.  Figure 
63  illustrates  the  corners  of  the  regressor  space  are  not  covered.  It  is  not  just  that 
there  are  not  data  points  there.  At  10  deg  AoA  and  25  deg/s  pitch  rate,  there  are  not 
data  points  either,  but  there  is  not  a  huge  inaccuracy  in  the  comparison  maneuvers  at 
that  value.  The  predictions  are  still  accurate  because  of  the  data  points  that  bound 
it.  It  would  make  sense  then  in  order  to  ensure  the  entire  regressor  space  is  covered, 
the  training  maneuver  should  encompass  the  regressor  space.  While  it  seems  to  be 
common  practice  in  the  literature  [15]  to  constrain  the  training  maneuver  within 
the  same  regressor  space  as  the  model  prediction  space,  there  is  no  reason  why  the 
training  maneuver  can  not  be  extended  beyond  that. 


Figure  63.  Edge  problem  shown  in  TM4  regressor  space 
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4.7  Lessons  Learned  Maneuver 


In  consideration  of  the  results  from  Section  4.6.1,  a  decidedly  different  approach 
was  taken  for  the  final  maneuver.  Most  of  the  previous  research  has  pointed  towards 
completing  training  maneuvers  within  specified  regressor  space  bounds.  At  the  very 
least,  the  literature  [17]  [15]  makes  the  claim  to  stay  within  the  ranges  of  the  training 
maneuver,  which  could  be  taken  to  mean  that  staying  within  max/min  Q  and  AoA 
would  be  sufficient.  Unfortunately,  it  is  entirely  possible,  as  this  research  can  show, 
that  it  is  possible  to  fill  the  regressor  space  and  still  overextend  the  models  within  the 
chosen  regressor  space.  By  extending  the  training  maneuver  outside  of  the  regressor 
space  ‘box’,  there  is  a  much  better  opportunity  to  ensure  model  accuracy  in  the 
entirety  of  the  regressor  space 

The  addition  of  a  single  sine  oscillation  to  the  beginning  of  the  maneuver  will  allow 
for  a  guaranteed  inclusion  of  the  entire  regressor  space  by  the  training  maneuver.  A 
single  sine  oscillation  at  constant  amplitude  and  frequency  forms  a  circle  on  the  Q 
vs  AoA  regressor  map.  Hitting  the  four  corners  of  the  regressor  space  ensures  there 
will  be  data  points  constraining  the  fit  in  the  entire  regressor  space.  Then,  the  “best 
practice”  DC  chirp  maneuver  is  used  to  fill  the  space  to  desirable  levels.  Figures  64 
and  65  illustrate  TM6. 

4.7.1  Lessons  Learned  Maneuver  Results 

The  RSP  values  for  TM6  versus  the  other  maneuvers  are  shown  in  Table  20.  While 
TM6  does  not  return  the  best  values  when  compared  against  the  other  maneuvers,  it 
does  have  values  in  the  same  range  as  the  other  ‘good’  maneuvers.  As  was  discovered 
with  TM5,  it  is  hard  to  individually  manipulate  a  single  RSP.  A  best  effort  has  to 
be  made  to  observe  changes  with  changing  inputs.  The  results  from  the  sixth 
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Figure  64.  Training  maneuver  6:  AoA  and  Q  vs  time 
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Pitch  Rate  (Q) 


Table  20.  RSP  values  for  all  TMs 


RSP 

Description 

Max  or  Min 

TM1 

TM2 

TM3 

TM4 

TM5 

TM6 

RSPn 

LQ  Q  std 

Min 

0.2047 

0.0024 

1.1168 

0.0158 

0.0234 

0.1842 

rsp17 

Time 

Max 

0.9143 

0.8851 

0.2293 

1.0005 

1.2660 

0.9670 

RSPU 

Kurtosis  HQ  Q 

Min 

24.574 

5.7243 

120.33 

21.999 

12.396 

37.5211 

rsp7 

OQ  std 

Min 

0.0178 

0.0057 

0.0488 

0.0160 

0.0116 

0.0128 

training  maneuver  are  shown  in  Table  21,  which  shows  the  raw  TM6  results  color 
coded  ’podium  style.’ 


Table  21.  TM6  model  metrics 


COM 

R2 

CL 

NRMSD 

R2 

cD 

NRMSD 

R2 

CM 

NRMSD 

COMla 

0.9823 

0.0334 

0.9748 

0.0430 

0.8716 

0.0755 

COMlb 

0.9721 

0.0419 

0.9624 

0.0562 

0.8198 

0.1328 

COMlc 

0.9880 

0.0247 

0.9812 

0.0351 

0.8817 

0.0890 

COM2 

0.9740 

0.0371 

0.9729 

0.0484 

0.8722 

0.1038 

COM3 

0.9911 

0.0244 

0.9735 

0.0284 

0.9074 

0.0568 

COM4 

0.9328 

0.0525 

0.9871 

0.0377 

0.9349 

0.0699 

COM5 

0.9594 

0.0466 

0.9682 

0.0590 

0.8250 

0.1452 

From  the  table,  it  can  be  seen  that  TM6  did  very  well  at  the  high  Q  values. 
However,  TM6  seemed  to  struggle  with  the  low  Q  maneuvers.  The  reason  for  this 
can  be  seen  in  Figure  66.  The  fit  to  the  data  is  good  in  the  high  Q  range;  however, 
the  fit  drops  off,  by  the  end  of  the  maneuver  which  is  in  the  middle  of  the  regressor 
space. 

The  Kestrel  MVP  in  Figure  67  was  constructed  using  the  same  SIDPAC  software 
as  used  in  this  research.  However,  due  to  the  organized  fashion  of  the  coefficient  results 
when  an  entire  aircraft  is  taken  into  consideration,  compared  to  a  1  ft  chord  airfoil,  the 
constructed  model  is  much  more  accurate  to  the  training  data.  The  difference  is  likely 
the  major  contributor  to  the  error  seen  by  TM6.  The  choice  of  the  desired  regressor 
space  may  have  been  too  ambitious  for  an  airfoil  from  the  start  of  the  current  research. 
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At  the  very  large  AoAs  and  Qs  generated  in  the  sixth  maneuver,  the  difference  in  flow 
phenomena  could  be  too  large  for  the  model  to  be  able  to  fit  the  entire  maneuver. 
It  would  be  interesting  to  contrast  the  use  of  a  RBF  (radial  basis  function)  model  in 
this  situation.  RBF  models  have  the  characteristic  of  exactly  matching  the  training 
data.  However,  the  theory  behind  the  generation  of  TM6  should  still  stand,  if  applied 
to  a  more  applicable  problem. 


Figure  66.  Fit  problem  shown  in  TM4  regressor  space 

Despite  TM6  not  performing  as  well  as  hoped,  there  is  a  lesson  learned.  Referenc¬ 
ing  Figure  68,  it  can  clearly  be  seen  the  edge  of  the  regressor  space  error  is  removed. 
TM6  is  very  accurate  in  that  region,  where  the  other  models  are  inaccurate.  How¬ 
ever,  TM6  has  problems  predicting  the  low  Q,  low  AoA  region  at  the  beginning  of 
the  comparison  maneuver. 
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Figure  67.  Fit  for  an  F-16  training  maneuver  completed  by  AFSEO  [11] 


Figure  68.  Model  results  COM3  Cd  with  TM6 
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V.  Conclusions  and  Future  Work 


5.1  Overview 

Identification  of  stability  issues  early  during  the  design  process  is  an  immediate 
need  for  the  Air  Force.  The  inability  to  detect  these  issues  has  led  to  costly  overruns 
and  late  design  changes  for  several  aircrafts.  Late  design  changes  lead  to  ad  hoc  fixes 
that  are  not  ideal  for  the  aircraft  function.  CFD  has  been  proposed  as  a  tool  for 
engineers  to  get  high-fidelity  static  and  dynamic  results.  CFD  has  a  large  computa¬ 
tional  cost  associated  with  it,  and  for  CFD  to  identify  these  trouble  areas,  requires 
an  immense  amount  of  simulations  in  a  variety  of  conditions.  System  identification 
has  recently  been  introduced  as  a  method  to  bring  CFD  and  S&C  together  by  reduc¬ 
ing  the  computational  cost.  The  CFD  SID  process  hinges  upon  the  decision  of  an 
appropriate  training  maneuver.  The  results  of  this  research  allow  engineers  to  make 
smarter  decisions  in  the  training  maneuver  development. 

This  research  includes  a  variety  of  topics  including:  CFD,  system  identification, 
S&C  and  statistics.  The  methods  and  actions  required  to  solve  the  Navier-Stokes 
equations  with  turbulence  models  to  get  high-fidelity  solutions  was  discussed.  The 
need  for  accurate  stability  and  control  parameters  has  been  demonstrated.  Previous 
research  was  considered,  and  the  lessons  herein  were  presented  and  applied  to  the 
current  work. 

A  C-Grid  topology  was  applied  to  a  NACA  64A010  airfoil,  and  a  method  for 
determining  grid  and  time  step  independence  was  presented  and  completed.  The  grid 
was  determined  to  be  time  and  grid  density  independent,  and  the  resulting  grid  and 
timestep  were  selected  to  move  forward.  When  dynamic  maneuvers  were  simulated, 
the  grid  was  in  need  of  further  refinement  to  eliminate  stability  issues. 

Several  different  training  and  comparison  maneuvers  were  created  and  simulated 


109 


using  CREATE-AV’s  flow  solver  Kestrel  and  were  evaluated  using  regressor  space 
parameters.  Multivariate  polynomial  models  were  created  from  the  training  maneu¬ 
ver  data.  The  coefficient  predictions  made  by  these  models  were  compared  against 
coefficients  obtained  from  the  full  CFD  simulation  of  comparison  maneuvers  and  suit¬ 
able  metrics  were  calculated.  Those  metrics  and  the  regressor  space  parameters  were 
fed  into  the  statistical  analysis  program  JMP®and  correlation  values  were  calcu¬ 
lated.  Those  correlation  values  were  used  to  refine  the  regressor  space  parameters 
and  shape  an  additional  training  maneuver.  Lessons  learned  from  conducting  the 
above  evaluation  led  to  the  creation  of  a  six  training  maneuver. 

5.2  Conclusions 

The  goal  of  this  research  was  to  investigate  methods  by  which  to  measure  training 
maneuvers  before  the  full  CFD  computational  simulation  and  to  link  the  results  of 
this  evaluation  process  to  future  performance.  Fulfilling  this  goal  would  create  a 
greater  understanding  of  how  to  devise  and  create  training  maneuvers.  While  the 
results  from  this  research  do  not  point  to  a  definitive  rules  for  creating  a  accurate 
training  maneuver,  the  results  are  suggestive  of  techniques.  Using  visual  methods, 
there  is  a  noticeable  difference  between  training  maneuvers  with  better  high  Q  RSPs 
and  high  Q  performance,  as  was  the  case  with  TM1  and  TM4  for  Cl  in  the  high  Q 
comparison  maneuver  COMlc.  For  the  low  Q  maneuver,  this  link  was  demonstrated 
with  TM1  providing  the  best  results  for  Cd- 

Correlation  values  between  the  metrics  and  RSPs  were  calculated  to  determine 
which  RSPs  were  most  predictive  of  favorable  TM  performance.  The  analysis  showed 
that  the  way  in  which  the  points  were  distributed,  as  measured  by  standard  deviation, 
was  important.  The  results  of  this  analysis  led  to  the  creation  of  a  total  of  36  RSPs. 
A  large  number  of  those  set  up  to  measure  how  the  points  are  distributed  in  the 
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regressor  space.  RSPn ,  RSPn,  RSP3e  and  RSP7  were  determined  to  provide  insight 
into  future  maneuver  performance.  RSPn  is  a  measure  of  the  rate  at  which  the 
training  maneuver  covers  unique  area  of  the  regressor  space.  RSPn  is  a  measure  of 
the  low  pitch  rate  standard  deviation  in  pitch  rate.  RSP^q  is  a  measure  of  skewness  in 
high  pitch  rate  in  pitch  rate.  RSP7  is  a  measure  of  zero  pitch  rate  standard  deviation. 

Using  the  selected  regressor  space  parameters,  a  fifth  training  maneuver  (TM5) 
was  created  with  specific  RSPs  in  mind.  The  creation  of  this  maneuver  proved  that 
it  is  hard  to  modify  all  these  RSPs  via  one  maneuver.  For  many  of  the  RSPs  it 
is  abstract  to  think  how  modifying  the  training  maneuver  motion  will  result  in  the 
desired  changes  in  the  RSPs.  TM5  was  generated  such  that  the  selected  RSPs  indi¬ 
cated  favorable  results.  The  results  TM5  were  favorable,  with  TM5  producing  the 
best  comparison  metrics  for  several  comparison  maneuvers  and  coefficients;  however, 
it  did  not  return  the  best  metrics  for  every  coefficient  and  every  maneuver.  A  type 
of  error  due  to  the  edge  of  the  regressor  space  was  observed  even  for  TM5. 

Lessons  learned  over  the  course  of  this  research  provided  the  needed  insight  for  the 
creation  of  TM6.  Rather  than  restrict  the  training  maneuver  to  the  confines  of  the 
defined  regressor  space,  providing  data  to  bound  the  extremes  of  the  box  was  decided 
to  be  an  avenue  worth  pursuing.  TM6  did  not  return  the  best  overall  model;  however, 
it  does  shed  light  on  the  limitations  of  MVP  models.  The  choice  of  a  regressor  space 
for  this  project  may  have  been  overly  ambitious,  with  much  of  the  regressor  space 
being  highly  unsteady.  Extending  that  range  of  AoA  and  Q  experienced  by  the  airfoil 
even  further  may  have  resulted  in  the  poor  fit  to  the  training  data  shown  in  Chapter 
IV. 

From  this  research  it  has  been  shown  that  RSPs  are  important.  Using  selected 
RSPs  as  a  guide,  an  additional  maneuver  (TM5)  was  created  that  was  able  to  out¬ 
perform  maneuvers  generated  without  the  RSPs  in  a  variety  of  comparison  scenarios. 
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This  research  has  shown  it  is  important  how  the  the  data  is  distributed  in  the  regres¬ 
sor  space  through  RSPs  that  measure  distribution  of  data  as  opposed  to  quantity  of 
data.  This  research  has  also  shown  that  pitch  rates,  such  as  Q  are  as  important,  if 
not  more  than  the  angles  such  as  angle  of  attack  through  the  RSPs. 

5.3  Future  Work 

While  much  ground  has  been  covered  in  this  work,  there  are  always  more  avenues 
that  could  be  explored.  The  first  possible  avenue  is  to  apply  the  higher  order  statistics 
using  the  previous  system  of  calculating  distribution  RSPs,  using  cells  instead  of 
number  of  hits.  While  the  author  tried  to  use  a  variety  of  methods  and  statistics 
for  quantifying  the  training  maneuvers,  it  was  not  an  exhaustive  effort.  Continued 
research  into  methods  and  statistics  could  be  fruitful.  Further,  utilizing  RBF  models 
particularly  in  regions  where  MVP  models  have  trouble  fitting  the  training  data,  as 
in  TM6,  could  produce  more  accurate  results. 

A  further  area  of  research  could  be  implementing  a  optimization  routine  to  the 
RSPs.  The  first  step  would  be  to  identify  a  cost  function.  Time  would  be  an  ap¬ 
propriate  choice.  Then  the  RSPs  would  need  to  be  set  up  as  constraints.  The  best 
choice  would  likely  to  be  determine  an  ‘acceptable’  level  for  the  RSPs  and  set  append 
the  RSPs  as  inequality  constraints.  Additional  constraints  would  be  angle  and  rate 
constrictions.  Second  derivative  constraints  may  be  necessary  to  ensure  a  continu¬ 
ous  path.  The  challenging  piece  would  be  including  abstract  concepts  like  standard 
deviation  and  skewness  into  the  mix. 

The  next  step  then  is  to  progress  into  3  dimensions.  The  lessons  learned  from  this 
research  should  be  applied  for  the  lateral  directional  coefficients  as  well.  Further,  the 
results  need  to  be  extended  for  other  flow  conditions.  This  research  was  run  at  Mach 
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0.5  at  10,000  ft.  What  may  be  important  for  model  prediction  accuracy  could  depend 
on  the  flow  conditions. 

Grid  generation  is  a  topic  that  is  brought  up  often  in  the  CFD  community  and  is 
acknowledged  by  practitioners  to  be  exceedingly  challenging.  There  are  methods  in 
the  literature  for  producing  a  static  grid;  however,  prescribed  motion  problems  require 
a  different  process.  While  the  technology  for  computing  more  advanced  problems 
improves,  methods  for  generating  grids  for  those  problems  needs  to  keep  pace. 

The  newest  version  of  Kestrel,  Version  3.1.1,  offers  the  ability  to  use  control  sur¬ 
faces  in  the  flow  solution.  This  offers  an  interesting  opportunity  to  further  integrate 
S&C  and  CFD.  While  it  is  possible  to  gain  insight  into  stability  and  control  using 
non-deflecting  grids  and  problems,  it  would  be  an  immense  value  to  the  S&C  field 
to  incorporate  derivatives  such  as  Cmge  which  is  change  in  pitch  rate  due  to  elevator 
deflection.  Furthering  CFD  and  SID  interaction  to  include  moving  control  surfaces  is 
a  very  exciting  prospect  for  future  research. 

This  research  has  shown  that  there  is  merit  to  quantifying  the  pertinent  charac¬ 
teristics  of  potential  training  maneuvers  in  order  to  allow  for  better  informed  training 
maneuver  generation.  While  adopting  a  best  practice  maneuver  may  produce  ade¬ 
quate  results,  by  researching  and  defining  what  characteristics  are  needed  for  a  train¬ 
ing  maneuver,  this  work  represents  a  step  forward  toward  introducing  a  standardized 
process. 
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TRAINING  MANEUVER  EVALUATION  FOR  REDUCED  ORDER  MODELING 
OF  STABILITY  &  CONTROL  PROPERTIES  USING  COMPUTATIONAL 

FLUID  DYNAMICS 


A.  Matlab@Pre  and  Post  processing  Files 


A.l  Spectral  Analysis  M-file 

1  */,2LT  Craig  C  Porter 

'/.Frequency  Analysis  MATLAB  File  in  Support  of  ROM  Thesis  Work 
'/.Derived  from  M-file  provided  by  Dr.  Keith  Bergeron 

clear 
6  clc 
close 

7. .  .  . 

7, Factors  that  need  to  be  changed 

11 

file_namel  =  ’ CGRID3T1 . coef f 1 ; 

cd  ../Data/St  70Change  depending  on  location  of  data 
datal  =  dlmread ( f ile.name 1  ,  ’  ’  ,20,0)  ; 

Iteration=datal (: ,1) ; 

16  Time=datal(: ,2) ; 

A0A  =  datal ( :  ,3)  ; 


BETA=datal ( : ,4) 

CAXIAL 

=datal ( 

,5) 

CN0RMAL 

=  dat al  ( 

,6) 

21 

CLIFT 

=  datal  ( 

,7) 

CDRAG 

=  datal  ( 

,8) 

CSIDE 

=  datal  ( 

,9) 

CPITCH 

=  datal  ( 

,  10) 

CR0LL 

=  datal  ( 

,11) 

26 

CYAW  =  datal ( :  ,  12)  ; 

N=max ( Iteration) ; 

cutoff  =  2000;  '/.if  desired  removes  startup  iterations 

del_time  =  Time  ( N ) -Time  (N -  1 )  ;  '/.assumes  constant  delT 

70L2  =  2048;  70  Length  of  signal  (truncated  to  a  ... 

power  of  2,  i.e.,  2048)  THIS  IS  WHAT  YOU  NEED  TO  CHANGE 

31 

70Default  is  CN0RMAL  for  FFT  ,  change  below  if  neccesary 
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%  Beware  line  below  is  correct,  i.e.,  L2  =  a  power  of  2! 

i =1 ; 

36  while  (N- cutof f ) >2~ i 
i=i+l ; 
end 

L2  =2  ~ ( i - 1 )  ; 

41  for  i=  l:N-cutoff 

amplitude  1 (i)=CN0RMAL (i  +  cutoff )  ; 

end 


amp.averagel  =  mean ( ampl itude 1 ) ; 

46 

for  i  =  l:N-cutoff 

y2(i)  =  ampl itude 1 ( i )- amp.aver age  1 ; 

end 


51 

for  i 


end 


1 : N- cutoff 
time  ( i ) 


( i  —  1 ) *del_t ime  ; 


56  f igur e  (  1 ) 

plot (time,  y2,  ’ color  ’,  [0  1  1]  ,’ Marker square LineWidth 

,.5)  ' 

hold  on 

xlabel(’Time  (s)’) 
ylabel ( ’ Amplitude  ’  ) 

61  t itle (’ Normal  Force’) 


1 1 1 1 1 1 1 7. 1 1 1 1 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 
7,  FFT 

7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 


66  Fs2  =  l/del_time;  7»  Sampling  frequency 

NFFT2  =  2~  nextpow2  ( L2 )  ;  7«  Next  power  of  2  from 

length  of  y 
n2  =  NFFT2 ; 

yfft2  =  fft (y2 , NFFT2) /L2 ; 

f  2  =  Fs2/2*linspace (0 , 1 , NFFT2/2  +  1)  ; 

71  power2  =  2* abs ( yf f t 2 ) ; 


7.  Plot  single-sided  amplitude  spectrum, 
f igure (2) 

plot  ( f 2  ,2*abs (yfft2 (1 : NFFT2/2+1) ) ) 

76  title (’ Single -Sided  Amplitude  Spectrum  of  Normal  Force’) 

xlabel (’ Frequency  (Hz)’) 
ylabel (’ Magnitude  of  Normal  Force’) 
axis  (  [0  500  0  1] ) 
axis  ’ auto  y  ’ 
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81 


l 


86 


91 


96 


101 


Identify  dominant  frequency 
temp2 ( 1 )  =  0.0; 

temp2  (2)  =  0.0; 

temp_ j  2  =  1; 
temp_power2  =  power2 ; 
max_pwr2 ( 1 : n2 /2)  =  0.0; 

max_f r eq2 ( 1 : n2 /2)  =  0.0; 

for  i  =  1 : n2 /2 

for  j  =  1 : n2 /2 

if  (temp_power2 ( j )  >  temp2(l)) 

temp2(l)  =  t emp_power 2 ( j ) ; 
temp2 (2)  =  f  2 ( j )  ; 

t emp_ j 2  =  j  ; 

end  ; 

max_pwr2(i)  =  temp2(l); 

max_freq2(i)  =  temp2(2); 
end  ; 

temp2 ( 1 )  =  0.0; 

temp2  (2)  =  0.0; 

temp_power2 ( t  emp_ j 2 ) =0 . 0 ; 

end 


'/.Print  top  6  frequency’s  to  screen 
106  for  i  =  1 :  6 ; 

i,  Amplitude  =  max_pwr2(i),  Frequency  =  max_freq2(i) 
end 

di sp  (  1  ’  ) 

di sp ( 1  Aver age  Ampltidue’) 

111  disp ( amp_averagel  ) 
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A. 2  Model  Generation  Script 


7.  Craig  Porter  TM  Model  Creation  Scipt 
7.7.  Load  Data 
3  clc 

clear  all 
close  all 

File= ’ TM6  1  ; 

8  file_namel  =  strcat(File,  ’.coeff  ’); 

70file_name2  =  strcat(File,  ’.flow’); 
file_name3  =  strcat(File,  ’.motion’); 

cd  /home / af it enl /gae 13m/ cport er /Thes i s /TM/ TM6b /  “/Change  depending 
on  location  of  data 

13  7.  “/.If  using  v2 

7.  datal  =  dl mr e ad (file_nam el,  ’’,5519,0); 

7.  7. data2  =  dlmread(file_name2,’’, 5519,0); 

7.  data3  =  dlmre  ad  (  f  i  le_name3  ,  ’  ’  ,  5519,0); 

18  "/.If  using  3.1.1 

datal  =  dlmread ( f ile_name 1  ,  ’  ’  ,5520,0)  ; 
data3  =  dlmread(file_n am e3,’’, 5522,0); 

TIME=datal ( : ,2) ; 

23 

“/.For  3.1.1  Runs,  refernce  frame  screwup  Q  is  actually  -R 
Q=-data3 ( : , 13) ; 

7,  7, For  v2  runs 
28  7.  Q  =  data3  (  :  ,  12)  ; 

A0A=datal ( : ,3) ; 

CLIFT  =  datal ( :  ,7)  ; 

CDRAG=datal ( ; ,8) ; 

33  CPITCH  =datal ( :  ,10)  ; 

ITERATI0N  =  datal ( :  ,  1)  ; 
dt  =  1 . 25  e -5 ; 

QDQT=deriv(Q,dt) ; 

A0AD0T  =  deriv(A0A,dt)  ; 

38 

clear  datal  data3 

7.7.  Create  Model 

43  close  all 
clc 

x  (  :  ,  1 )  =  ADA  ;  “/.ALPHA 
x  (  :  ,  2)  =Q  ;  7.Q 
48  x ( : , 3) =QDQT ; 
z  =  CLIFT  ; 
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NORD  =5 ; 

nord= [NORD  NDRD  NORD]; 
maxord  =5 ; 

53  sig2=0; 
auto=l ; 

Iplot  =1  ; 
i var  =0 ; 
bvar  =0 ; 

58  maxopt=0; 

[y,p,  ip  ,  crb  ,pse  ,xp  ,a,  ia,psi]=mof  (x,z  ,  nord  , maxord  ,  sig2  ,  auto  ,  Iplot  . 

ivar , bvar , maxopt ) ; 

[xr,xlab]=polygen(x,ip) ; 

63  f igure (2) 

plot (TIME , z , TIME , y) 
xlabel(’Time  (s)’) 
ylabel ( ’ Coef f i cent  of  Lift’) 
legend ( ’ CFD  ’  ,  ’ Model ’ ) 
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A. 3  Initial  RSP  M-file[5] 


function  [metric  metrica  metricq  as  qs  z]=Metric... 

2  (pos2,rate2, bounds .steps .metricsteps , maxNO , pmO) 

'/.extract  AoA  and  Q  bounds  of  regressor  space 
aextrema  =  bounds (1  ,  : )  ; qextrema  =  bounds (2 ,  : )  ; 

'/.extract  step  sizes  for  AoA  and  Q  for  regressor  space  ... 
discretization 

7  ast ep  =  st eps  (  1 )  ; qst ep  =  st eps (2)  ; 

7,  Debugging  variables 

ametstep=metricsteps (1) ; qmetstep=metricsteps (2) ; 

12  '/.All  AoA  values 

arange=aextrema(l) :astep:aextrema(2) ; 

7.  All  Q  values 

qrange=qextrema(l) :qstep:qextrema(2) ; 

17 

7, Create  2D  Mesh  for  discretized  regressor  space 
[as,qs]=meshgrid(. . . 

mean ( arange (1 :2) ) : astep :mean(arange (length(arange) -1: length ( 
arange ))),... 

mean ( qrange (1 :2) ) : qstep :mean (qrange ( length (qrange) -1: length ( 
qrange ) ) ) ) ; 

22 

°/0Z-axis  of  3D  plot  of  regressor  map  (number  of  data  points  per  . 
cell ) 

z=zeros (size (as) ) ; 

'/.Make  AoA  ’  s  (pos  variables)  and  Q  values  same  length 
27  pos3=zer os ( 1 , length ( pos2 )- 1 ) ; 
for  i =2 : length ( pos2 ) 

pos3(i-l)=mean(pos2(i-l: i)) ; 

end 

32  '/.Loop  through  AoA  1  s  to  mark  which  cell  each  point  lies 
for  i = 1 : length ( rat e2 ) 

if  rate2 (i) >qextrema (1)  &&  rate2 ( i ) <qextrema (2) ; 

'/.Find  which  AoA  cell 
for  j = 1 : length ( arange ) -1 
37  apos  =  j  ; 

if  pos3 ( i )<= arange ( j + 1 ) 
break 

end 

end 

42  '/.Find  which  Q  cell 

for  j = 1 : length ( qrange ) -1 
qpo s  =  j  ; 

if  r at e2 ( i ) <=qr ange ( j + 1 ) 
break 
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47 


end 


end 

’/.Increment  "number  of  data  points  per  cell"  counter 
z ( qpos , apos)=z(qpos , apos ) + 1 ; 

end 

52  end 

asteps=size (z ,2) / ametstep ; 

'/.Examine  the  evenness  of  the  data  in  AoA 
counta=zeros ( asteps ,1) ; 

57  totala=zeros ( asteps ,1) ; 

totala (1) =size (z , 1) * ametstep ; 
for  i  =  l : size  (z  ,  1) 
k  =  l ; 

for  j =1 : size (z  ,  2) 

62  if  z  (  i  ,  j  )  ~  =  0 

count a (k)=counta(k)+l; 

end 

if  mod ( j , amet st ep ) ==0 

totala (k) =totala  ( 1)  ; 

67  k=k+l; 

end 

end 

end 

metrica=counta . / totala ; 

72 

'/.Examine  the  evenness  of  the  data  in  Q 
qsteps=size (z , 1) / qmetstep ; 
countq=zeros ( qsteps ,1) ; 
totalq  =  zeros (qsteps  ,1)  ; 

77  totalq (1) =size (z , 2) *qmetstep ; 
for  i=l:size(z,2) 
k  =  l; 

for  j  =1 : size (z  ,  1) 
if  z  (  j  ,  i  )  ~  =  0 

82  count q  (k)  =  count q  (k)  + 1  ; 

end 

if  mod ( j , qmet st ep ) ==0 

totalq (k) =totalq  ( 1)  ; 
k=k+l ; 

87  end 

end 

end 

metricq=countq . / totalq ; 

°/o  ***  **  Metric  1  -  total  percentage  of  cells  with  data  points 
92  metric (1) =100*  sum (countq) / sum (totalq)  ; 

7.  Ex  amine  the  data  points  on  the  boundaries 
count  =0 ; 
total =0 ; 

97  j  =  size (z  ,  2)  ; 

for  i  =  l : size  (z  ,  1) 
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if  z  (i  ,  1)  ~=0 

count  =  count  + 1 ; 

end 

102  total=total+l ; 

if  z  (i  ,  j  )  ~=0 

count  =  count  + 1 ; 

end 

total=total+l ; 

107  end 

j  =  size (z  ,  1)  ; 
for  i=2 : size (z , 2) -1 
if  z  (1  ,  i)  ~=0 

count  =  count  + 1 ; 

112  end 

total=total+l ; 
if  z  (  j  ,  i)  ~=0 

count  =  count  + 1 ; 

end 

117  total=total+l ; 

end 

'/.Find  the  location  of  zero  pitch  rate 
for  i = 1 : length ( qrange ) -1 
122  zeroqpos=i; 

if  0 <= qrange ( i + 1 ) 
break 

end 

end 

127  ”/.*****  Metric  2  -  total  percentage  of  boundary  cells  with  data  .. 
point  s  *  *  *  *  * 

metric (2) =100*  count /total ; 

°/o*****  Metrics  3  &  4  -  Evenness  of  data  in  AoA  and  Q  respectively 
metric (3 : 4) = [std (metrica/max (metrica) )  std (metricq/max (metricq) ) ] 

132 

'/.Isolate  the  data  points  at  0  Q 
zsteps=size (z ,2) / qmetstep ; 
countz  =  zeros (zsteps  ,1)  ; 
k  =  l; 

137  i=zeroqpos; 

for  j=l:size(z,2) 
if  z  (i  ,  j  )  ~=0 

countz (k)=countz (k)+z (i , j ) ; 

end 

142  if  mod ( j , qmet st ep ) ==0 

k=k+l ; 

end 

end 

'/.Cap  max  number  of  data  points  in  a  given  cell  (if  desired) 

147  countzl2=0; 

for  i = 1 : length ( count z ) 
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if  countz ( i ) >maxN0 

countz (i)=maxN0 ; 

end 

152  if  countz  (i)  >0 

countz 12=countz 12+1 ; 

end 

end 

V, *****  Metric  5  -  total  percentage  of  0  Q  cells  with  data  points 

157  metric (5) =100*  count z 1 2/length ( countz )  ; 

“/,  *  *  *  *  *  Metric  6  -  Avg  #  data  points  per  cell  at  0  Q  ***** 
metric (6)=me an (countz) ; 

162  “/„*****  Metric  7  -  evenness  of  data  at  0  Q  ***** 
metric (7)=std( countz /max ( countz ) ) ; 

°/0Now  isolate  region  of  "low"  Q  (excluding  region  of  0Q) 
countz22  =0 ; 

167  for  mi = 1 : length ( pmO ) 

count z2 =[z( zero qp os -pmO (mi)  :zeroqpos-l  ,  :)  ;z(zeroqpos  +  l:  ... 
zeroqpos+pmO(mi)  ,:)]  ; 

‘/.Cap  max  number  of  data  points  per  cell  (if  desired) 
for  i  =  1 : s ize ( countz2  ,  1 ) 

for  j =1 :  size ( countz2  ,  2) 

172  if  countz2 ( i , j ) >maxN0 

countz2 (i , j ) =maxN0 ; 

end 

if  countz2 ( i , j ) >0 

count z22=countz22+l ; 

177  end 

end 

end 

“/,  *  *  *  *  *  Metric  8  -  total  percentage  of  "low"  Q  cells  with  data  ... 
point  s  ***** 

metric (6  +  2*mi) =100*  count z 22 / (size ( count z 2  ,l)*size( countz 2  ,  2) ) 

182 

°/o*****  Metric  9  -  Avg  #  data  points  per  cell  at  "low"  Q  ***** 
metric (7+2*mi)=mean(mean(countz2) ) ; 

'/.There  may  be  NaN’s  in  metric  calculation  after  normalization 
-  change 
187  '/.these  to  0 

for  i  =  1 : s ize ( countz2  ,  1 ) 

met8(i)  =  std(countz2(i,  :)/max( countz 2 ( i  ,  : ) ) )  ; 

end 

j  =i ; 

192  for  i=l : length (met8) 

if  i snan ( met 8 ( i ) ) 
met82 ( i ) =0 ; 

else 

met82(i)=met8(i) ; 
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197 


end 


end 

°/t  *  *  *  *  *  Metric  10  -  evenness  of  data  at  "low"  Q  ***** 
metric (8+2*mi)=mean(met82) ; 

end 

202 

7,  Isolate  region  of  "high"  Q 
mil =9+2*mi ; 

countz3  = [z (1 : zero qp os -pmO (mi)  ,:)  ;z(zeroqpos  +  pmO(mi)  : size (z  ,  1)  ,  :  )  ] 
countz32  =0 ; 

207  for  i  =  1 : s ize ( count z3  ,  1 ) 

for  j =1 : size ( countz3  ,  2) 
if  countz3 ( i , j ) >0 

countz32=countz32+l ; 

end 

212  end 

end 

°/o*****  Metric  11  -  total  percentage  of  "high"  Q  cells  with  data  . 
point  s  *  *  *  * 

metric (mil) =100*  count z32 / (size ( count z 3  ,l)*size( count z 3  ,  2) )  ; 

217  °/0*****  Metric  12  -  Avg  #  data  points  per  cell  at  "high"  Q  ***** 
metric (mll+1) =mean (mean ( count z 3 ) ) ; 

'/.There  may  be  NaN  ’  s  in  metric  calculation  after  normalization  -  . 

change 
“/» these  to  0 

222 

for  i  =  1 : s ize ( count z3  ,  1 ) 

metl0(i)=std(countz3(i,  :)/max(countz3(i,  :)))  ; 

end 

j=i; 

227  for  i  =  1 : length ( met  1 0 ) 
if  i snan (met  10 ( i ) ) 
met  102 ( i ) =0 ; 

else 

metl02(i)=metl0(i) ; 

232  end 

end 

°/o*****  Metric  13  -  evenness  of  data  at  "high"  Q  ***** 
metric (mil +2) =mean (met  102 )  ; 

237  end 
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A. 4  Final  RSP  calculation  M-file 


function  [rsp , as ,qs ,z]=RSP(pos2 ,rate2 ,  that ,dt , bounds , steps ,maxN0 
pmO  ) 

3  “/extract  AoA  and  Q  bounds  of  regressor  space 
aextrema  =  bounds (1  ,  : )  ; 
qextrema  =  bounds  (2,  :) ; 

“/.extract  step  sizes  for  AoA  and  Q  for  regressor  space  ... 
discretization 
8  ast ep  =  st eps  (  1 )  ; 
qstep  =  steps  (2)  ; 

“/.All  AoA  values 

arange=aextrema(l) :astep:aextrema(2) ; 

13 

“/„  All  Q  values 

qrange=qextrema(l) :qstep:qextrema(2) ; 

“/.Create  2D  Mesh  for  discretized  regressor  space 
18  [as , qs] =meshgrid ( . . . 

mean ( arange (1 :2) ) : astep :mean(arange (length(arange) -1: length ( 
arange ))),... 

mean (qrange (1 :2) ) : qstep :nean (qrange ( length (qrange) -1: length ( 
qrange ) ) ) ) ; 

“/.Z-axis  of  3D  plot  of  regressor  map  (number  of  data  points  per  . 
cell ) 

23  z=zeros ( size (as) ) ; 

NOBS  =0 ; 

“/.Make  AoA  ’  s  (pos  variables)  and  Q  values  same  length 
pos3=zeros (1 , length (pos2) -1) ; 

28 

for  i =2 : length ( pos2 ) 

pos3(i-l)=mean(pos2(i-l: i)) ; 

end 

33  “/.Loop  through  AoA’s  to  mark  which  cell  each  point  lies 
for  i =1 : length ( rate2 ) 

if  r  at  e2  (  i  )  >=  qext  r  ema  (  1 )  kk  r  at  e2  (  i  )  <=  qext  r  ema  (2 )  ; 

if  pos3 ( i ) >=aextrema  ( 1)  kk  pos3 ( i ) <= aextrema (2)  ; 

38  “/.Find  which  AoA  cell  it  is  in 

for  j = 1 : length ( arange ) -1 

if  pos3 (i) <=arange ( j +1) 
apos  =  j  ; 
break 

43  end 

end 

“/.Find  which  Q  cell 
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48 

53 

58 

63 

68 

73 

78 

83 

88 

93 


for  j =1 : length ( qrange ) -1 

if  r at e2 ( i ) <=qr ange ( j + 1 ) 
qpo  s  =  j ; 
break 

end 

end 

"/Increment  "number  of  data  points  per  cell"  counter 
z(qpos ,apos)=z(qpos , apos ) + 1 ; 

N0BS=N0BS+1 ; 
end 

end 

end 

asteps  =  size (z  ,  2)  ; 
qsteps  =  size (z  ,  1)  ; 

"/Examine  the  evenness  of  the  data  in  AoA  and  Q  and  count 

r spa=zeros ( asteps ,1) ; 

r spq=zeros ( qsteps ,1) ; 

counth=0 ; 

for  i=l : asteps 

for  j  =1 : qsteps 

"/For  every  loop  sum  the  number  of  hits  for  each  AoA 
rspa(i)=rspa(i)+z(j ,i) ; 

"/Every  loop  sum  the  number  of  hits  for  each  Q 
rspq(j)=rspq(j)+z(j ,i) ; 

“/Keeping  track  of  how  many  cells  are  hit 
if  z ( j  ,  i  )  >0 

counth=counth+l ; 

end 

end 

end 

"/Examine  the  data  points  on  the  boundaries 
countb  =0 ; 
total =0 ; 
j  =asteps ; 

"/Left  and  right  boundaries 
for  i=l : qsteps 

if  z  (  i  ,  1)  ”  =  0 

countb =countb+l ; 

end 

t  ot  al  =  t ot  al  + 1 ; 
if  z  (i  ,  j  )  ~=0 

countb=countb+l ; 

end 

t  ot  al  =  t ot  al  + 1 ; 

end 

j  =qsteps ; 
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98 

103 

108 

113 

118 

123 

128 

133 

138 

143 


%Top  and  bottom  boundaries ,  making  sure  to  not  double  count  the  . 
corners 

for  i=2:asteps-l 
if  z  (1  ,  i)  ~  =  0 

countb=countb+l ; 

end 

total=total+l ; 
if  z  (  j  ,  i)  ~  =  0 

countb=countb+l ; 

end 

t  ot  al  =  t ot  al  + 1 ; 

end 

"/.Find  the  location  of  zero  pitch  rate 
for  i=l : qsteps 

if  0 <= qr ange ( i  + 1 ) 
zeroqpos  =  i  ; 
break 

end 

end 

'/.Count  the  total  number  of  data  points  at  0  Q 
countz=zeros ( asteps ,1) ; 
j  =  zeroqpo  s ; 
for  i=l : asteps 

countz(i)=z(j  ,i)  ; 

end 

'/.Cap  max  number  of  data  points  in  a  given  cell  (if  desired) 
countzl2=0 ; 
for  i=l : asteps 

if  countz ( i ) >maxN0 

countz (i)=maxN0 ; 

end 

if  countz ( i ) >0 

countz 12=countz 12+1 ; 

end 

end 

'/.Isolate  region  of  "low"  Q  (excluding  region  of  0Q) 
countlq=0 ; 

countlq2= [z (zero qp os -pmO : zero qp os -1 , : ) ; z (zeroqpos+1 : zeroqpos+ 
pmO  ,  :  )  ]  ; 

"/.Cap  max  number  of  data  points  per  cell  (if  desired) 
for  i  =  l : size ( countlq2  ,  2) 

for  j = 1 :  size ( count lq2  ,  1 ) 

if  countlq2 ( j , i) >maxN0 

countlq2 (j , i) =maxN0 ; 

end 

if  countlq2 ( j , i) >0 

count lq= count lq+1 ; 

end 
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148 


end 


end 

s izelq  =  s ize ( count lq2  ,  1 )* s ize ( count lq2  ,2)  ; 

'/Isolate  region  of  "high"  Q 

counthq2= [z (1 : zero qp os -pmO  ,  :)  ; z(zeroqpos+pmO : size (z  ,  1)  ,  :  )  ]  ; 

153  counthq=0; 

for  i  =  1 : s ize ( counthq2  ,  2) 

for  j =1 : s ize ( counthq2  ,  1 ) 

if  counthq2 ( j , i ) >maxN0 

counthq2 (j , i)=maxN0 ; 

158  end 

if  counthq2 ( j  ,  i )  >0 

counthq=counthq+l ; 

end 

end 

163  end 

sizehq=size ( counthq2  ,l)*size(counthq2  ,2)  ; 

%RSP  calcs 

rspalq  =  zeros ( size ( countlq2  ,2)  ,1)  ; 

168  rspqlq=zeros(size( count lq2 , 1 ) , 1 ) ; 
rspahq  =  zeros (size ( counthq2  ,2)  ,1)  ; 
rspqhq  =  zeros (size ( counthq2  ,  1 )  ,  1 )  ; 
for  j =1 : size ( countlq2  ,  1) 

for  i  =  l : size ( countlq2  ,  2) 

173  rspalq(i)=rspalq(i)+countlq2(j,i); 

'/.Every  loop  is  summing  the  number  of  hits  for  each  Q 
rspqlq(j)=rspqlq(j) +countlq2 ( j , i) ; 

end 

end 

178  for  j =1 : size ( counthq2  ,  1) 

for  i  =  l : size ( counthq2  ,  2) 

rspahq(i)=rspahq(i)  +  counthq2(j  ,i)  ; 

'/.Every  loop  is  summing  the  number  of  hits  for  each  Q 
rspqhq(j)=rspqhq(j)+counthq2(j ,i) ; 

183  end 

end 

'/.Pearson  Cal  c 

rho 1 =zer os ( length ( pos3 )  ,1)  ; 

188  rho2=zeros (length (pos3) ,1) ; 
rho3  =  zer os ( length ( pos3 )  ,1)  ; 
mean3=mean (pos3) ; 
mean2=mean (rate2) ; 
for  i = 1 : length ( pos3 ) 

193  rhol(i)=(pos3(i)-mean3)*(rate2(i)-mean2) ; 

rho2(i)=(pos3(i)-mean3) ~2  ; 
rho3(i)=(rate2(i)-mean2) ~ 2 ; 

end 

198  °/t  *  *  *  *  *  RSP  1  -  total  percentage  of  cells  with  data  points  ***** 
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rsp (1) =100*counth/ (asteps*qsteps) ; 

°/o  *  *  *  *  *  RSP  2  -  total  percentage  of  boundary  cells  with  data  points 
rsp  (2) =100*countb/total ; 

203 

*  *  *  *  *  RSP  3  &  4  -  Evenness  of  data  in  AoA  and  Q  respectively  ... 
rsp(3:4)=[std(rspa/(asteps))  std(rspq/(qsteps))] ; 

°/o*****  rsp  5  -  total  percentage  of  0  Q  cells  with  data  points  ... 

208  rsp (5) =100* countzl2/asteps ; 

“/,  *  *  *  *  *  RSP  6  -  Avg  #  data  points  per  cell  at  0  Q  ***** 
rsp (6) =mean ( countz )  ; 

213  “/„*****  RSP  7  -  evenness  of  data  at  0  Q  ***** 
rsp (7) =std (countz/asteps) ; 

°/0*****  rsp  8  -  total  percentage  of  "low"  Q  cells  with  data  points 
rsp (8) =100*  count lq/ (sizelq)  ; 

218 

°/o*****  rsp  9  -  Avg  #  data  points  per  cell  at  "low"  Q  ***** 
rsp (9) =mean (mean ( countlq2 ) ) ; 

°/o*****  rsp  10  &  11  -  evenness  of  data  at  "low"  Q  ***** 

223  rsp (10: 11) =[std(rspalq/ size ( count lq2 ,2))  std(rspqlq/size( count lq2 . 

,D)]; 

°/o*****  rsp  12  -  total  percentage  of  "high"  Q  cells  with  data  ... 
point  s  ***** 

rsp (12) =100*  counthq/ (sizehq)  ; 

228  °/0*****  RSP  13  -  Avg  #  data  points  per  cell  at  "high"  Q  ***** 
rsp (13) =mean (mean ( counthq 2 ) ) ; 

°/o*****  rsp  14  &  15  -  evenness  of  data  at  "high"  Q  ***** 
rsp(14:15)=[std(rspahq/size(counthq2 ,2))  std(rspqhq/size( counthq2 . 

,D)]; 

233 

°/o*****  rsp  16  -  Training  Maneuver  Effectiveness  ***** 
r sp ( 16 ) =N0BS / length ( rat e2 ) ; 

°/o*****  rsp  17  -  Time  Ef  f  e  ct  i  vene  s  s  1  ***** 

238  r sp ( 17 ) =r sp ( 1 ) /max ( that ) ; 

°/o*****  rsp  18  -  Time  Effectiveness2  ***** 
rsp(18)=rsp(l)*dt /max (that) ; 

243  °/„*****  RSP  19  &  20  &  21  Time  Eff  0Q  ,  LQ  ,  HQ  **** 
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rsp (19:21)=[rsp(5)/max (that)  rsp (8) /max (that)  rsp (12) /max (that)] 

°/t*****  RSP  22  &  23  -  Kurtosis  Alpha  k  Q  ***** 
rsp (22 : 23) = [kurtosis (rspa)  kurtosis (rspq) ] ; 

248 

°/o*****  rsp  24  &  25  -  Skewness  Alpha  &  Q  ***** 
rsp (24 : 25) = [skewness (rspa)  skewness (rspq) ] ; 

°/o*****  rsp  26  -  Pearson  product -moment  correlation  coefficient 
253  rsp  (26 )  =  sum  (rhol)/sqrt  (  sum  (  rho2  )  *  sum  (  rho3  )  )  ; 

°/o*****  rsp  27  k  28  Kurtosis  and  Skewness  OQ  ***** 
rsp (27 : 28) = [kurtosis ( countz )  skewness ( countz ) ] ; 

258  “/„*****  RSP  29  k  30  Kurtosis  LQ  AoA  /  Q  **** 

rsp (29 : 30) = [kurtosis (rspalq)  kurtosis (rspqlq)] ; 

%*****  rsp  31  &  32  Skewness  LQ  AoA  /  Q  **** 
rsp (31 : 32) = [skewness (rspalq)  skewness (rspqlq)] ; 

263 

°/o*****  rsp  33  &  34  Kurtosis  HQ  AoA  /  Q  **** 
rsp (33 : 34) = [kurtosis (rspahq)  kurtosis (rspqhq)] ; 

°/o*****  rsp  35  &  36  Skewness  HQ  AoA  /  Q  **** 

268  rsp (35 : 36) = [skewness (rspahq)  skewness (rspqhq) ] ; 


end 
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B.  Figures  with  TM5  included 


Figure  69.  Model  Results  COMla  Cl 
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Figure  70.  Model  Results  COMla  Cd 


Figure  71.  Model  Results  COMla  Cm 


131 


Figure  72.  Model  Results  COMlb  Cl 


Figure  73.  Model  Results  COMlb  Cd 
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Figure  74.  Model  Results  COMlb  Cm 


Figure  75.  Model  Results  COMlc  Cl 
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Figure  76.  Model  Results  COMlc  Cd 


Figure  77.  Model  Results  COMlc  Cm 
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Figure  78.  Model  Results  COM2  Cl 


Figure  79.  Model  Results  COM2  Co 
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Figure  80.  Model  Results  COM2  Cm 
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Figure  81.  Model  Results  COM3  Cl 


137 


Figure  82.  Model  Results  COM3  Co 


Figure  83.  Model  Results  COM3  Cm 
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Figure  84.  Model  Results  COM4  Cl 


Figure  85.  Model  Results  COM4  Co 
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Figure  86.  Model  Results  COM4  Cm 


Figure  87.  Model  Results  COM5  Cl 
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Figure  88.  Model  Results  COM5  Co 


Figure  89.  Model  Results  COM5  Cm 
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